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We consider the effects of new light species on the Cosmic Microwave Background. In 
the massless limit, these effects can be parameterized in terms of a single number, the 
relativistic degrees of freedom. We perform a thorough survey of natural, minimal models 
containing new light species and numerically calculate the precise contribution of each of 
these models to this number in the framework of effective field theory. After reviewing the 
relevant details of early universe thermodynamics, we provide a map between the parameters 
of any particular theory and the predicted effective number of degrees of freedom. We then 
use this map to interpret the recent results from the Cosmic Microwave Background survey 
done by the Planck satellite. Using this data, we present new constraints on the parameter 
space of several models containing new light species. Future measurements of the Cosmic 
Microwave Background can be used with this map to further constrain the parameter space 
of all such models. 



Contents 



I. Introduction 

II. Methodology 

A. Relativistic Species and the CMB 

1. Silk Damping 

2. Early Integrated Sachs- Wolfe Effect 

B. Early Universe Thermodynamics 

C. Decoupling, Recoupling, and the Redistribution of Entropy 

D. Relevant Decoupling Temperatures 

E. Renormalization Group Effects 

F. Big Bang Nucleosynthesis 

G. Summary 

III. Models 

A. Spin-0: Goldstone Boson 

B. Spin-i; Light Fermion 
f. Gauge Interactions 

2. Dipole Moments 

3. Four- Fermion Interactions 

4. Sterile Neutrinos 

5. Axinos 

C. Spin-f : Gauge Boson 

1. Kinetic Mixing and Gauge Interactions 

2. Dipole Moments 

D. Spin- 1 : Gravitino 

E. Spin-2: Graviton 

F. Models Without New Light Species 

IV. Conclusions 

A. Notations and conventions 



2 



B. Details of code 



66 



Acknowledgments 



References 



72| 

72 



I. INTRODUCTION 

The Cosmic Microwave Background (CMB) is one of the only probes we have of physics 
in the early universe. Through a detailed mapping of anisotropics in the temperature of 
those photons which decoupled from visible matter in the era of recombination, we are able 
to determine the relativistic energy density in that era. From this, we gain information 
about the number of light species in our universe. In the massless limit, we can accomplish 
this by a fit to only one number, the relativistic degrees of freedom jl-5]. This parameter 
is often expressed in terms of an effective number of neutrinos, N e ff, defined such that 
in the Standard Model (SM) of particle physics N e ff is roughly the number of neutrino 
generations. Beyond the Standard Model (BSM) physics models which contain new light 
species with masses (9(eV) or less can contribute to this measurement. Consequently, we 
have new terrain in which to test the SM through its prediction of g* = 3.38, corresponding 
to an N eff of 3.046 |g-lfj]. 

There has been a statistically insignificant but consistent excess in the measured value 



of g* 



17H20j. Prior to the results from the Planck satellite, the most precise reported 



measurement was g* = 3.69 ± 0.16, corresponding to N e ff = 3.71 ± 0.35 19|, coming from 
a combination of data from the South Pole Telescope (SPT) and the Wilkinson Microwave 
Anisotropy Probe (WMAP). A similar excess is present in measurements from the Atacama 
Cosmology Telescope (ACT) 18) . Very recently, however, the Planck collaboration released 
the first results from its measurement of CMB anisotropies, obtaining a result of g* = 
3.50 ± 0.12, corresponding to N e ff = 3.30 ± 0.27 [2l|. Future Planck results will continue 
to improve the precision of this measurement, with a projected final g* sensitivity of ±0.09 



22 



. l23l| . In addition, future measurements of the polarization of the CMB are pro jected to 
constrain g* to within ±0.02, corresponding to constraints on N e ff of ±0.044 [23| . We are 
entering an era of being able to contrast the SM prediction for g* with the predictions of 
BSM physics models containing new light species to an unprecedented precision. 
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The power of this probe of new physics is that in any BSM theory containing new species 
with masses < 0.1 eV which were once in thermal equilibrium with the SM, the effect of 
these species is contained in a single number, the correction Ag* to the SM prediction for 
<?*. Therefore, a map from the parameters of a BSM model to the number A g* can be 
constructed in order to determine the consistency of regions of the parameter space with 



the measured value of g*. Although useful approximations of such a map exist [24], [25J, we 
are entering the exciting era of precision cosmology experiments, and consequently it has 
become imperative to form precise theoretical predictions. The subject of this paper is the 
precise numerical computation of this map of model parameters to Ag* for a wide variety of 
natural, minimal BSM theories containing new light or massless species. We approach this 
problem in a largely model-independent effective field theory framework to fully characterize 
the effects of all such models. 

Although there are other existing constraints on new light species pre sent in the early 
universe coming from the study of Big Bang Nucleosynthesis (BBN) 26M28| , this probe does 
not have the same resolving power as the Planck satellite. Unlike BBN, Planck and future 
polarization experiments have the power to probe the actual values of the couplings of new 
light species to the SM, as we shall demonstrate in this work. Even in the absence of a signal 
for new physics from future experiments, the results of this work provide new constraints on 
the couplings of SM species to new light particles which are competitive with, and sometimes 
even surpass, existing constraints from other areas of physics. This establishes a new arena 
for testing the predictions of BSM physics models with new light species. 

The recent results for g* from Planck are in tension with independent measurements of 



the Hubble expansion rate today 



2l|. Specifically, combining those measurements with the 



results from Planck leads to a preference for higher values of g* than quoted above. Therefore, 
these results are not capable of confirming or rejecting the hypothesis of new light degrees 
of freedom being present in the early universe. Regardless, in order to demonstrate the 
constraining power of measurements of the CMB, we proceed as if this tension were not 
present. Motivated by the Planck results given above, we proceed by supposing that values 
of g* > 3.74 (N e ff > 3.84) are excluded at the 95% confidence level. We interpret our results 
in this framework in order to illustrate how further data could be utilized. 

In section [Til of this paper, we review the relevant details of the determination of g* 
using the CMB, as well as details of thermodynamics in an expanding universe, providing a 
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framework for the rest of the paper. In section HTT1 of this paper, we discuss all BSM physics 
models compatible with our criteria of naturalness and minimality. Specifically, we discuss 
the parameters which provide the interaction strength between various fields in the SM and 
the new light species present in the model. We present the current experimental constraints 
on each of these scenarios, as well as our findings for the contribution of each new light 
species to g* as a function of the parameters in the underlying theory. We also interpret 
the viable parameter space of each model in terms our aforementioned interpretation of the 
recent results from the Planck satellite, placing additional constraints on theories using this 
new CMB data. 

II. METHODOLOGY 

We study the effects of adding new light or massless particles to the SM on the evolution 
of the universe and the CMB. Specifically, we investigate new particles which at some time 
in the early universe were in equilibrium with the SM and decouple prior to recombination. 
Translating between additional fields in the Lagrangian and the measurement of the effective 
number of light degrees of freedom, g*, requires a detailed analysis of the quasi-thermal 
evolution of the universe. The effects of new light degrees of freedom depend on both when 
and how they decouple from the thermal bath. As we shall see, a direct measurement of 
anisotropies in the CMB then leads to a resultant measure of g* at recombination. 

In this section, we first review how light species predominantly affect the CMB, namely 
via Silk damping and the early integrated Sachs- Wolfe (ISW) effect. We also review the 
thermodynamics of the early universe, as well as the effects of decoupling and other non- 
equilibrium events. We then discuss the range of decoupling temperatures which can signif- 
icantly impact the CMB. Finally, we briefly review the most important existing constraint 
on new light degrees of freedom, namely their effect on Big Bang Nucleosynthesis. 

A. Relativistic Species and the CMB 

The early universe was not perfectly homogeneous, but instead had small perturbations 
in the distribution of energy density, which are currently believed to be seeded by inflation. 
These regions of under- or overdensity correspond to small perturbations in the metric away 
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from the pure Friedmann-Robertson- Walker form. CMB anisotropies provide a direct mea- 
surement of these early universe perturbations, whose distribution and structure are sensitive 
to the thermodynamic conditions leading up to recombination. The CMB therefore gives us 
insight on the properties and structure of the universe in its infancy. The measurement of 
g* using the CMB is performed through a precise determination of the expansion rate, H, in 
the era of recombination. The relationship between H and g* arises because the expansion 
rate is determined solely by the total energy density, p, and the curvature. Increasing the 
value of g* at a fixed temperature leads to a larger overall p, which then leads to more 
rapid expansion. Silk damping is sensitive to the value of H leading up to and during 
recombination, while the early ISW effect is affected by the evolution of H once photons 
are effectively free-streaming, which lasts from recombination onwards. For more detailed 



and thorough explanations of these effects than those presented here, consult [29l-l3l| and 
references therein. 



1. Silk Damping 

Prior to recombination, protons, electrons, and photons interacted very strongly to form 
a tightly-coupled plasma. Despite the high frequency of interactions, the mean free path 
for photons was nonzero, and photons were able to diffuse outward. The rate of photon 
diffusion grew as the protons and electrons combined into hydrogen, up until the point of 
last scattering. The overall diffusion scale at the end of recombination is therefore pre- 
dominantly determined by the mean free path during recombination and the duration of 
recombination. The diffusion of photons results in a partial thermalization of the baryon- 
photon plasma, damping any inhomogeneities on scales smaller than the photon diffusion 



e in turn leads to a damp- 
32], above some multipole 



length. This reduction of inhomogeneities below some length sea 
ing of temperature anisotropies, commonly called Silk damping 
moment la- A larger value for H then leads to a decrease in the amount of time available 
for this diffusion, restricting the damping to smaller angular scales and reducing the mag- 
nitude of the damping. An increase in would therefore lead to reduced Silk damping, or 
equivalently a larger damping moment. 

Any map between the predicted diffusion length and the precise value for Id is sensitive to 
experimental uncertainty in the angular distance to the last scattering surface. In practice, 
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it is simpler to remove this uncertainty by considering the ratio of Id to the smaller sound 
horizon moment l s . This sound horizon arises independently of photon diffusion, due to the 
spread of inhomogeneities in the baryon-photon plasma. These oscillations propagate at the 
corresponding speed of sound, setting an acoustic oscillation length scale at recombination. 
The addition of new light species reduces the time for these inhomogeneities to spread, 
which increases the value of l s , in addition to the increase in Id- These two processes, photon 
diffusion and sound wave propagation, have different time dependencies. This difference 
results in an increase of the ratio l s /ld as H grows, leading to damping of more of the 
acoustic peaks, despite the fact that the overall damping has been reduced. 



2. Early Integrated Sachs-Wolfe Effect 

Following recombination, photons propagate freely without scattering but pass through 
points of matter over- or underdensity. If the gravitational potential of these inhomogeneities 
is constant in time, there is no net effect on the CMB photons. However, if the gravitational 
potential has any time-dependence, the photons will experience some net loss or gain in 
energy as they pass through a single gravitational perturbation and be red- or blueshifted as 
a result. The alteration of CMB anisotropies due to time-dependent gravitational potentials 



is the ISW effect [33]. 



The evolution of gravitational potentials is determined by the expansion rate, which 
depends on the overall particle content. In a universe consisting solely of nonrelativistic, 
pressureless matter, the competing effects of gravitational clustering and universe expansion 
cancel, such that potentials are time- independent. However, any nonnegligible pressure 
alters the expansion rate such that the potentials do evolve with time. There are therefore 
two points in time at which the ISW effect could contribute to the CMB. The first occurs 
when the universe contains a nonnegligible radiation density, which is the case immediately 
following recombination. This alteration to the CMB shortly after its formation is commonly 
referred to as the early ISW effect. The second era corresponds to the point at which the 
vacuum energy becomes a significant fraction of the total energy density. This second case, 
which begins near modern times, is the late ISW effect. 

Unsurprisingly, new light species increase the radiation energy density following recombi- 
nation, altering H and enhancing the early ISW effect. Specifically, the presence of additional 
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FIG. 1: Projected CMB anisotropy power spectrum for three different values of g* (or equivalently 
N e ff). The addition of new light degrees of freedom increases the height of the first peak through 
the early ISW effect and decreases the height of later peaks through Silk damping. The power 
spectrum, and therefore these effects, are measured by multiple observational experiments, such as 
the Planck satellite. These spectra were calculated using CAMB [34j, l35j. The magenta, blue, and 
orange curves (dark gray, black, and light gray curves, when viewed in black and white) correspond 
to an N e jf of 3, 4, and 5, respectively. 



species causes gravitational potentials to evolve more rapidly, resulting in more substantial 
red- and blueshifts to CMB photons passing through these evolving potentials. On very 
small scales, photons will pass through multiple such potentials, and the net effect cancels. 
However, the potentials rapidly become time-independent, such that photons are unable to 
pass through multiple large-scale perturbations before this effect ends. An increase to g* 
therefore enhances the variance in temperature anisotropies on angular scales corresponding 
to the largest structures immediately following recombination. The size of the largest struc- 
tures at this point coincides with the acoustic horizon, such that the early ISW effect leads 
to an increase in the first acoustic peaks of the CMB. In practice, this effect is measured by 
comparing the height of the first acoustic peak to that of latter peaks. 

The effects of Silk damping and early ISW, which can be seen in Figure HJ are comple- 
mentary means of measuring H, and therefore g*, near recombination. However, they are 
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still sensitive to two different points in time. Silk damping probes H prior to and during 
recombination, while the early ISW effect is sensitive to H immediately after recombina- 
tion. This has two important consequences for constraints on light species. The first is that 
experiments which focus on precision measurement at smaller values of I, such as WMAP, 
are sensitive mainly to the early ISW effect. The resulting constraints are therefore more 
limited experimentally by the effects of cosmic variance. Experiments which instead focus 
on anisotropics at larger I, such as ACT, SPT, and Planck, are predominantly sensitive to 
Silk damping and are less affected by cosmic variance. 

The second consequence to note is that particles with masses near the temperature scale 
of recombination (~ 0.1 — 1 eV) will potentially contribute very different signals to these two 
sets of experiments. The detection of such species would involve a detailed analysis of each 
individual effect, rather than a simple fit to all the experimental data. While we consider 
massless particles for the majority of this work, we will return to this possibility for the case 
of sterile neutrinos. 



B. Early Universe Thermodynamics 

As mentioned earlier, CMB measurements of light species are predominantly sensitive to 
the relativistic energy density, which is characterized by an effective number of relativistic 
degrees of freedom g*. For more details on the material discussed in this subsection, see 



36M38| . We can then define g* in terms of p re i, the energy density of all relativistic species, 



and a reference temperature T, which we take to be the photon temperature (T = T 7 ), 

Prel = 9*7^T 4 . (1) 

The light species content of the SM, which consists of photons and neutrinos, can then 
be used to make a prediction for the measured value of g* at recombination. This prediction 
can be written in the form 

7 /T\ 4 7 /4\ 4/3 

9* = g, + ^ N eff{f) = 2 + 8' 2 ' Neff {u) ' (2) 

where g 7 = 2 is the number of spin degrees of freedom for photons and g v = 2 is the 
number of spin degrees of freedom for a single neutrino species. The factor of | is due 
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to the effect of Fermi-Dirac statistics on energy density, and T v is the calculated neutrino 
temperature assuming neutrinos instantaneously decouple from the rest of the SM at T ~ 
MeV. The parameter N e ff is the effective number of neutrino species. This historically 
defined parameter, which is 3.046 1 for the SM, is often used to parametrize the effect of any 
light species other than photons on g*. The contribution of neutrinos and any new light 2 
species to g* is given solely by N e ff. Any measured deviation from the SM prediction of 
= 3.38 would then indicate the need for new physics. 

This paper calculates the full contribution Ag* of new light species present in a large 
number of beyond the SM theories. This contribution to the relativistic degrees of freedom 
is found by calculating the energy density of new species near the point of recombination. 
The contribution can also be expressed as a change to N e ff as 



To find the energy density of a light species at recombination, we must track the evolution 
of its phase space density f{t,p). This form for the distribution function relies on the 
assumption that the universe is homogeneous and isotropic 3 . This phase space density is 
defined such that the particle number density is given by 



where g is the number of spin degrees of freedom for that particular species. We must first 
determine / at high temperatures, when the new species is in equilibrium with the SM, 
then calculate the changes to / as the universe expands and cools, with various species 

1 This effective number of neutrinos is denned such that if neutrinos truly did decouple instantaneously, 
N e ff would be 3. However, detailed calculations have shown this to not be the case, and the actual energy 
density of neutrinos is slightly larger than in the instantaneous decoupling approximation due to their 
interactions with annihilating electrons. This then results in the slightly larger predicted value for N e ff. 
For details on these calculations, see 



2 By light, we mean m <§C eV. The contribution of species with masses ~ eV is more complicated, as we 
shall discuss later. 

As discussed earlier, the universe is in fact not perfectly homogeneous or isotropic, and the distribu- 
tion functions therefore have some spatial and directional dependence. However, these deviations are 
quite small in magnitude, and any resulting correction to the CMB is below the experimental resolution. 
Consequently, any inhomogeneities and anisotropics in the distribution functions are negligible for our 
purposes. 




(3) 




(4) 
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annihilating or decoupling. 

As the universe expands, the evolution of each individual phase space density is controlled 
by both the rate of expansion H = -, where a(t) is the scale factor for the expanding universe, 
and the rate of interaction with the other particle species. This dependence is expressed 
using the Boltzmann equation 

where p — \p\ and the collision functional C[f] accounts for changes to / due to interactions. 
If we assume that the dominant interactions will consist of 2-to-2 scattering, then C[f] for 
some new species X is defined as the sum over all such possible interactions involving X. If 
each interaction process is time-reversal invariant, 



c lfx\ = \ E / ( II 9s tSMe) ^^-Po^siMfnUxJUjJk), (6) 

X,i^j,k \s=i,j,k ^ s / 

with the squared amplitude \A4\ 2 averaged over the spins of both incoming and outgoing 
particles. The term S corresponds to a symmetry factor whose value is \ when j and k 
are identical particles (e.g. e~ and e~ in the final state, but not e~ and e + ), to avoid 
overcounting of states in the phase space integral, and is 1 otherwise. The function 
is the phase space weighting term 

n(/*, h f s , f k ) = fMi ± fx)(i ± h) - fxfi(i ± mi ± /*), (7) 

where the ± term is + for bosons (Bose enhancement) and — for fermions (Pauli exclusion). 
The collision terms therefore couple together the Boltzmann equations for various particle 
species. 

Let us first consider the simple case of C[f] = 0, where species are freely evolving due 
solely to the expansion of the universe. Changing variables to write the phase space density 
for each species as f(t, q), where q = ap, we can rewrite the Boltzmann equation as 
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Put in this simpler form, we can easily see that any function of only q satisfies the homo- 
geneous Boltzmann equation. The phase space density only changes due to the redshifting 
of p, with no additional time- dependence. This result implies that for a freely evolving 
species, the naive guess of a Fermi-Dirac or Bose-Einstein distribution for / breaks down 
as the temperature falls to scales comparable to the particle mass m, as those distributions 
explicitly depend on \/p 2 + m 2 and can no longer be written solely as functions of q. In the 
massless limit, however, these distributions are possible solutions, provided the temperature 
scales as T ~ a -1 . 

We can now consider the effects of interactions, such that C[f] is potentially nonzero, 
but assume a negligible particle asymmetry, such that any chemical potential // « and the 
phase space density for particles is equal to that of antiparticles 4 . The general distribution 
for a particle species can be expressed as 

^(*>P) = e E/T(t)+<p(t,p) ± I ■ ( 9 ) 

The function <f)(t,p) simply parametrizes any deviations from the standard Bose-Einstein 
or Fermi-Dirac distributions. Such deviations can arise in two distinct situations. 

The first occurs if any species has a nonzero mass m. Once T < m, inelastic collisions 
which produce the massive species will be heavily suppressed. However, elastic collisions 
can still maintain kinetic equilibrium between this massive species and the remainder of the 
SM. If we insert our general form for / into equation (J7|) for the case of elastic collisions, 
we see that Q({f}) is still zero, and equilibrium is still maintained, so long as <ft is solely a 
function of time and not momentum. The (f)(t) term then arises when any species becomes 
nonrelativistic but is still in kinetic equilibrium with the remaining particles. This correction 
is not isolated to the distribution function of the massive species, but 'infects' those of all 
other species. This effect occurs in the SM as T reaches the mass of the electron. In this 



4 For the case of a light boson in chemical equilibrium with the SM, it turns out this approximation is exact, 
provided the boson is its own antiparticle. A light fermion, however, could potentially have a nonzero 
H, but this would not generate any asymmetry in the SM sector, due to baryon and lepton number 
symmetries. The only observable consequence of a nonnegligible chemical potential would then be an 
increase to the energy density (and therefore Ag*) of our new species. This asymmetry would need to be 
generated in some UV completion, making this extension very model-dependent and beyond the scope of 
this current work. 
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case, electromagnetic interactions are strong enough such that both photons and electrons 



maintain pure Bose-Einstein and Fermi-Dirac distributions for the bulk o 



annihilation. However, a small but nonzero <fi(t) arises once T < ^ 



39, 



the net electron 



40J. However, these 



particular deviations are small and consequently irrelevant for our purposes. 

The second possible situation for such deviations occurs when a species begins to decouple 
from the remainder of the SM during the annihilation of some other species. In this case, 
the momentum-dependence of C[f] causes a full <ft(t,p) alteration to the distribution of the 
decoupling species. This effect arises solely when a species is leaving equilibrium and does 
not affect the distribution functions of the remaining species in the thermal bath, provided 
that the rate of self-interactions of the thermal bath is sufficiently large. In the SM, this 
situation occurs as neutrinos decouple from electrons. The annihilation of electrons alters 
the momentum dependence of the neutrino distribution, introducing a nonzero 0. In fact, it 
is precisely this term which causes the deviation of the SM N e ff prediction from 3 to 3.046. 

A much more detailed treatment of these effects, as well as the full evolution of species 
in equilibrium, can be found in J39]. For our purposes, the important fact is that standard 
thermal distributions will not fully suffice for any non-equilibrium behavior, specifically the 
decoupling or annihilation of a species. For these general phase space density must 

be numerically evolved in time to find the precise contribution to g* at lower temperatures. 
The focus of this work includes both decoupling and annihilation, necessitating our numerical 
treatment. 

So far we have treated the expansion of the universe as an independent process, but it is 
in fact coupled to the evolution of its particle content through the Einstein field equations. 
Assuming a flat, isotropic, and homogeneous universe, we obtain the Friedmann equations, 
which can be written in the form 



H> = —p, 

where p and P refer to the total energy density and pressure of the full particle content. 
The Boltzmann equations and Friedmann equations then combine to give a coupled set of 
integro-differential equations governing the full evolution of the early universe. 
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C. Decoupling, Recoupling, and the Redistribution of Entropy 

While a full solution to equations (JSJ) and ffTUj) is necessary to understand the detailed 
evolution of any species X and its exact contribution to gr*, we can first gain a qualitative 
understanding by considering the approximation of instantaneous decoupling. Once we have 
developed this conceptual framework, we will then turn to more precise statements about 
the complete evolution of distribution functions. 

In the instantaneous decoupling approximation, the point of decoupling can be found by 
comparing the rate of expansion H to the rate of interaction T x , defined as 

Tx= J2 !hpi( av } jk ^ Xh (11) 

j,k—>-X,i 

where (av) is the thermally-averaged cross-section for any interaction j, k — >■ X,i. This 
average cross-section can be formally defined as 



Mi,^,i=/( II ^ 7 ^^H27r)V(p m -p ou 05|^| 2 ^(l±/x)(l±/,). (12) 

J \s=X,i,j,k { ' s I UjUk 



Note that the symmetry factor S now includes an additional factor of | if the initial state 
consists of identical particles, as well as the original \ for an identical-particle final state. 
The full set of thermally- averaged cross-sections can be related to the collisional term C[f] 
via 



/ 9x , n ^ (2ir) 4 8 4 (pin-Pwt)C[f x }S = y~] {n j n k {av) j ^x,i-n x n i {av)x,i^j,k) ■ (13) 
j [Ztt) hj X £—r 

X,t,j,k 

Conceptually, V x corresponds to the rate of production per particle for species X. In the 
limit where all species are massless, this rate will scale as Tx ~ T n , where n is determined 
by the scaling dimension of the relevant interaction operators. As the universe expands, 
both H and Tx will decrease, though generically at different rates. If Tx decreases more 
quickly than H, then it is possible for a species originally in equilibrium to 'freeze out' and 
decouple from the remainder of the SM. 
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Conversely, if H decreases more quickly than Fx, a species originally out of equilibrium 
may actually recouple to the SM. In this case, however, X will generically not have the same 
temperature as the SM, if it even has a well-defined temperature, prior to recoupling. In- 
stead, the initial distribution will depend on any other particle content that could potentially 
couple to X, making this scenario very model-dependent. 

This point of instantaneous decoupling/recoupling is defined simply as the temperature 
at which Tx = H. It is common to assume that species are in full equilbrium prior to 
decoupling, then evolve freely immediately after freezing out. This approximate description 
is correct only when all relevant species are relativistic and originally in full equilibrium. 
However, if X decouples during other nonequilibrium processes, such as nonrelativistic an- 
nihilation, the full set of Boltzmann equations must be used. 

Once T drops below the mass of any particle, that species begins to annihilate away, with 
the number density quickly falling to a negligible amount. The entropy of the annihilating 
species is redistributed amongst the remaining interacting species, such that the temperature 
of all remaining species decreases less quickly than would be the case in free expansion. If 
X has decoupled from the SM prior to this annihilation, it will not participate in the 
resulting entropy redistribution, and therefore reaches a temperature lower than that of the 
SM following the annihilation. This behavior is almost the case for neutrinos, which begin 
to decouple from the SM at T ~ MeV, such that only those neutrinos with large momenta 
receive entropy from the electron annihilation. 

To determine the impact of these entropy redistributions, we need to track the relativistic 
entropy density s as a function of temperature. If the entropy density of all SM species in 
equilibrium (excluding X) was initially s when X instantaneously decoupled from the SM 
at temperature T , conservation of total entropy gives us the resulting temperature ratio 
following an entropy redistribution. This ratio can be expressed as a function of the entropy 
density s of all species in equilibrium at any future temperature T, 



In practice, because the entropy of annihilating species is only being distributed amongst 
relativistic species in full thermal equilibrium, it is much simpler and equivalent to instead 
use the relativistic degrees of freedom, rather than s/T 3 , to calculate the ratio 




(14) 
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Tx 
T 




.after \ 1 / 3 



(15) 



where g% e and gf ter are the relativistic degrees of freedom of all SM species in equilib- 
rium immediately before and after the entropy redistribution. This decrease in relative 
temperature also decreases the Ag* due to X, 



where Ag* is simply the initial contribution of X to g* at T . For multiple entropy redis- 
tributions, the overall ratio can be found simply by multiplying together the ratios from 
each individual redistribution, giving the full contribution of X to g*. 

This is precisely the method used to calculate the neutrino temperature in the instanta- 
neous approximation. Once neutrinos decouple, the only other species contributing to the 
entropy density are electrons and photons, and following the net annihilation of electrons, 
only photons remain. This method gives the ratio of temperatures 



If there were an additional species X which interacts strongly enough with the SM such 
that the entropy from electrons is evenly distributed amongst both the photons and X, then 
Tx = T after the electron annihilation, with no further redistributions to raise the photon 
temperature. The participation of species X in the electron entropy redistribution would 
affect our calculation for equation f lTTl) . due to the contribution of X to the entropy density. 
Since the value of T is already fixed by experiment, the presence of X increases the predicted 
neutrino temperature (relative to the photon temperature), resulting in 



In this particular case, the change in effective degrees of freedom due to the addition of 
X has two contributions. The first is the degrees of freedom in X (Ag*x), and the second is 
the degrees of freedom added by effectively warming the neutrinos. So for any species with 
Tx = T after the electron annihilation, 




(16) 




(17) 




(18) 
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Ag* = Ag* x + -j 

Again, this discussion has made the simplifying assumption of instantaneous decoupling. 
In general, we cannot simply use comparisons of Tx to H to determine the exact evolution of 
X if the species decouples during nonequilibrium processes. Our treatment must instead be 
made more precise by numerically solving the full set of coupled Boltzmann and Friedmann 
equations. The SM particles interact strongly enough to maintain equilibrium and can 
therefore be described by a single temperature T(t). The only exception to this is neutrinos, 
but at temperatures T > MeV, neutrinos also maintain equilibrium. Any new species X 
needs to be described by its phase space density fx{t,E). Finally, the expansion of the 
universe is described by the scale factor a(t). The evolution of these three functions can 
then be found by numerically solving the coupled set of three equations 



2Ag< 



x 



U + 2Ag* x J 



n 4/3 
11 J 



(19) 



H 2 = ^-(psm + P x), (20) 



3 

dpsM , dp x 
dt 



+ -%jr = -ZH{p S M + Psm + Px + Px), 



where the SM thermodynamic properties psu and Psm can be expressed solely in terms of 
T(t). For temperatures near the point of neutrino decoupling, the three neutrino species 
must then be described using three additional phase space densities, each with their own 
Boltzmann equation. 

For a given model of new light species X, the collisional term C[/x] can be found in 
terms of the model parameters, such as the coefficient A of the operator of nonrenormal- 
izable couplings in an effective theory. This interaction term then governs the process of 
decoupling X from the SM. Any SM annihilation and entropy redistribution that occurs 
after this decoupling reduces the change in effective degrees of freedom Ag* at the point 
of recombination. The contribution to g* for a specific model can be found by using the 
resulting fx near the point of recombination to calculate the energy density px- Solving this 
contribution in terms of generic couplings establishes a direct relationship between model 
parameters and Ag*. 
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D. Relevant Decoupling Temperatures 

The focus of our work is on those models which require detailed numerical analysis. We 
are interested in those species which were at some point in thermal equilibrium with the SM 
and will contribute significantly to g*. Specifically, we consider species which decouple during 
other SM events, such as the entropy redistribution resulting from the net annihilation of 
some massive species. 

One possibility is that a new light species decouples at temperatures above all relevant 
scales of the SM, specifically at T > TeV. At such high temperatures, all SM species would 
be present and relativistic. Any potential contribution of this new species at recombination 
would then be diluted by the large number of entropy redistributions as all SM species 
eventually annihilate away. Using equation f[T6"j) . combined with the total SM degrees of 
freedom g* = 106.75, we find that the largest possible contribution of a species with gx 
degrees of freedom which decouples above all SM energy scales is 



< 0Mg x . (21) 

Any additional species beyond those of the SM which redistribute their entropy after 
X decouples will only further suppress the contribution Ag*. These contributions are far 
below the experimental resolution on g*, so new light species must decouple at much lower 
temperatures to be detectable with the CMB. The full dependence of Ag* on the decoupling 
temperature for various particle types is shown in Figure [2j This functional dependence is 
calculated in the instantaneous decoupling approximation by using equation ffl6|) in combi- 
nation with of the SM as a function of temperature, which is shown in Figu re p 

As we see in Figure [2} for a species to be within the sensitivity of Planck 21], it must 
decouple at temperatures T < 200 MeV, which corresponds to the approximate scale of 
the QCD phase transition (see j^l] and references therein for details). Prior to this point, 
quarks and gluons are the relevant degrees of freedom for the QCD sector, such that the total 
number of SM degrees of freedom is g* — 61.75. As the universe cools to lower temperatures, 
the SM transitions to a regime where mesons and baryons are the appropriate degrees of 
freedom. Specifically, the relevant hadrons present below the QCD phase transition are 
pions and charged kaons, such that g* = 19.25. This significant reduction in the degrees 
of freedom results from the rapid annihilation or decay of any more massive hadrons which 



18 



bo 
<l 




Z 
< 



0.01 1 
Decoupling Temperature (GeV) 

FIG. 2: Additional light degrees of freedom Ag* at recombination for a new light species as a 
function of the decoupling temperature (in the instantaneous decoupling approximation), calculated 
using equation (|16j) . The contribution of various particle species is shown, specifically a real scalar 
boson (magenta), a Weyl fermion (blue), a real gauge boson (orange), and a Dirac fermion pair 
(green). The dashed line indicates the current sensitivity of the Planck observational experiment 
2l| . The gray region corresponds to the QCD phase transition, where the precise evolution of 
g*(T) for the SM is not well-understood. The provided values of Ag* should therefore only be 
interpreted qualitatively in that region. 



may have formed during the transition. The QCD phase transition therefore corresponds 
to a large redistribution of entropy into the remaining degrees of freedom, such that any 
species which decouples from the SM prior to the transition will not contribute significantly 
to the CMB. 

In principle, it is possible to discover species which decouple during the QCD phase 
transition, as those species could contribute values of Ag* above the experimental sensitivity. 

3ecause of, 



42 



44J and 



However, the precise details of this phase transition are not well-understood 
e.g., strong coupling effects, and this transition is an area of active study (see 
references therein). Consequently, we do not know how to make precise predictions for Ag* 
for species decoupling in this era. These computations are beyond the scope of our work, so 
we choose to restrict our focus to species which decouple after the QCD phase transition. 
For new species which do decouple immediately after this point, the calculation of Ag* is 
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FIG. 3: Effective degees of freedom in the SM as a function of temperature. The gray region cor- 
responds to the QCD phase transition, where the precise evolution of g*{T) is not well-understood. 
The provided values of g* should therefore only be interpreted qualitatively in that region. 



sensitive to whether the species couples to leptons or to quarks. Species which couple solely 
to leptons have a straightforward decoupling process, as all relevant interactions are suffi- 
ciently weakly renormalized. Species which couple to quarks will then couple to pions and 
kaons, whose couplings can be strongly renormalized. We must restrict ourselves to quark 
and meson couplings which involve conserved currents, as these are then protected against 
strong renormalization effects. For this set of couplings, we can still make precise predic- 
tions for the contribution of new light species which couple to quarks, even at temperatures 
immediately below the QCD phase transition. 



E. Renormalization Group Effects 

In this approach, we consider couplings between the SM and new light fields through 
higher-dimensional operators of the form A 4_d (9, where O is some operator comprised of 
both SM and new fields and d = dim(O). We consider the effective field theory of each 
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operator at very low energies (as low as ~ MeV). In order to match to the full UV theory 
which generates these operators, one should in principle treat A 4_d as a Wilson coefficient 
and run this coupling from the high energy theory down to the scale of interest using the 
renormalization group. We assume that this running has already been done when we write 
down our effective operators, such that we are working with the matched coefficient. In 
practice, this effect is not a large one, as we are always either in the weakly coupled regime 
or considering conserved current operators. In cases where the running begins at the GUT 
or Planck scale, however, large logarithms could in principle be generated which would need 
to be computed before our results can be utilized. 



F. Big Bang Nucleosynthesis 

Most models which include additional light degrees of freedom will have other model- 
dependent constraints, such as those from collider signals or various astrophysical observa- 
tions. Arguably the most important model-independent bound other than that of the CMB 
is that placed by Big Bang Nucleosynthesis (BBN). The measurement of the primordial 
relic abundance of light elements formed by BBN provides an independent probe of new 
light species, although at times earlier than recombination. While here we only give a brief 
summary of the relevant aspects of BBN, an excellent introduction to the topic can be found 



in [261428] . 



As the universe cooled to temperatures below the characteristic scale of nuclear binding 
energies (T < MeV), protons and neutrons began to bind together to form various light 
elements, such as deuterium, tritium, and helium. However, the large binding energy per 
nucleon of helium-4 ( 4 He) prevented any significant formation of further elements. The large 
majority of neutrons were therefore incorporated into 4 He. 

The abundances of the light elements, particularly 4 He, are therefore sensitive to the 
number density of neutrons at the start of BBN. When neutrons and protons were in full 
equilibrium, the number of neutrons relative to that of protons continued to fall due to their 
mass splitting. The neutron abundance is then determined by the point at which the weak 
interactions, which interconvert protons and neutrons, freeze out. A larger expansion rate 
results in earlier freezeout, which in turn leads to a larger number of neutrons and therefore 
more 4 He. 
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The precise value of H at the time of BBN, which would be increased by the presence of 
additional light species, therefore determines the relic abundance of 4 He. This abundance is 
often expressed in terms of the so-called 'helium mass fraction' 



Y P = 



4n 



He 



(22) 



n H + 4n H e 

Measurements of Yp therefore provide another means of constraining the relativistic en- 
ergy density of the early universe, though it is important to remember that these constraints 
apply at a different period of time than those placed by direct CMB measurements of g*. 
The SM prediction for the helium abundance is Yp = 0.2487 ± 0.0006 26|, and this predic- 



tion can be tested by both direct observation of helium and the detection of effects of Yp on 
CMB anisotropies. 

Multiple previous helium measurements have yielded results near Yp = 0.240 ± 0.006 



26l | . which are consistent with SM predictions, but two recent observational 



indicated a higher abundance of Y P = 0.2565 ± 0.0010 (stat) ±0.0050 (syst) 



studies have 



0.2561 ± 0.0108 



45[ and Y P 



461 ]. which are consistent with a larger rate of expansion. This in turn 



allows for the presence of new light species. In addition, combined CMB constraints from 
SPT and WMAP have measured Yp = 0.296 ±0.030 47J, and combined results from Planck 
and WMAP have measured Yp = 0.266 ± 0.021 [21]. These measurements are therefore 
currently incapable of either completely confirming or excluding the existence of new light 
species, but instead increase the importance of the precision CMB measurements of g* 
possible with future experiments. 

Lastly, it is important to note that there is tension between the SM prediction and ob- 
servational measurements of the abundance of lithium-7 ( 7 Li), with a lower observationally 
inferred primordial 7 Li abundance than that predicted by BBN. Unfortunately, this dis- 
crepancy is not immediately remedied simply by the presence of new light species, and the 
detailed model-building necessary to address this tension is beyond the scope of this pa- 



per. However, the 7 Li prob 



discovery of new physics 



48 



em does present another exciting opportunity for the possible 
452] - 



22 



G. Summary 

We have now introduced the framework necessary for the remainder of this paper. The 
focus of this work is the effects of light species in BSM theories on the CMB, which we 
determine by computing the energy density of the new light species at recombination. We 
specifically concern ourselves with species which were in thermal equilibrium with the SM 
and then decouple after the QCD phase transition, potentially during the annihilation of a 
SM species. Species which decouple from the SM before the QCD phase transition cannot 
be probed by the Planck satellite. The energy density of light species is calculated by numer- 
ically solving for the distribution functions of the new light species as the species decouple 
from the SM. Specifically, we numerically solve the coupled Boltzmann and Friedmann 
equations, found in equations §5§ and ffTUl) . in order to compute the potentially nonthermal 
distribution function given in equation ([9]). The distribution function immediately follow- 
ing decoupling can then be used to calculate the energy density at recombination, which 
determines g* using equation (JT|. 



III. MODELS 



In this section, we consider the set of models which can contribute to the CMB measure- 
ment of g* 5 . We mainly restrict ourselves to models where the additional degrees of freedom 
were in thermal equilibrium immediately following the QCD phase transition 6 . Such models 
must either contain new species with mass < eV or alter the neutrino energy density. While 
there are a very large number of possible models one could write down, we choose to restrict 
ourselves to those which are both minimal and natural. 

We consider a model to be minimal if it contains the smallest possible hidden sector in 
the low-energy theory. In particular, this restricts our discussion to models of elementary 
particles, ignoring the possibility of light composite states. We then direct our attention 



5 We only consider models with light degrees of freedom. It is possible to construct models where heav- 
ier species mimic the effects of light degrees of freedom through a nonzero presure resulting from non- 
equilibrium distribution functions |53| . 



54 



6 There are models where out-of-equilibrium effects such as decays generate a contribution to g* 
but no generic model-independent statements can be made about such scenarios, so we do not consider 
them in this work. 
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to the low-energy effective field theory (EFT) and ignore any additional particle content 
which may arise at higher energies, as these are irrelevant for our calculations. For each 
such model, we present a simple UV completion for the sake of motivation, but any UV 
theory in the same universality class will result in the same contribution to generalizing 
our results to a large number of possible theories. 

For this work, we define naturalness as technical naturalness. We therefore require that 
the size of quantum corrections not exceed the size of the physical observables in the theory, 
i.e. |^| < 1 for all parameters A, as large corrections require an artificial fine-tuning of 
parameters. 

A large number of potential models of light species are unnatural, due to large corrections 
to the mass of that new species. There are two predominant methods of suppressing quantum 
corrections to a particle's mass. The first method is the introduction of an additional 
symmetry which prohibits the existence of a mass term for that species. The second option 
is the use of strong dynamics in a hidden sector to generate large anomalous dimensions for 
mass terms, such that those terms become irrelevant operators, giving rise to a vanishing 
mass in the low-energy EFT. However, most models of the latter type tend to contain a 
relatively rich spectrum, violating our minimality principle. Although this is an interesting 
direction for future research, it is outside the class of models we consider. We therefore focus 
solely on theories of light species which contain a protective symmetry. 

The classes of possible new light species can be divided up by spin, as this restricts the 
protective symmetries available. We progress through each possible case, from spin-0 to 
spin-2, considering all minimal, natural models. For scalar fields, only a shift symmetry can 
preserve naturalness. Such fields are realized as Goldstone modes of a spontaneously broken 
global symmetry. Supersymmetry can also protect against loop-induced scalar masses, but 
only if supersymmetry is unbroken in the hidden sector, and even then, direct coupling 
to the SM as well as supergravity effects will communicate the breaking to the scalar. A 
sufficiently low supersymmetry breaking scale is in strong conflict with the spectrum of 
particles observed at colliders. Fermions can be protected by a chiral symmetry, and gauge 
bosons are protected by a gauge redundancy. We will employ these symmetries to generate 
minimal models with light degrees of freedom. Finally, theories with elementary spin-| and 
spin-2 particles are greatly restricted and are not good candidates for contributions to g*, 
as we shall review later in this section. 
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For each model, we then scan over all allowed couplings, numerically solving the Boltz- 
mann and Friedmann equations to calculate the full process of decoupling for any species 
which decouples during SM entropy redistributions. The details of our numerical approach 
can be found in Appendix [B] Using the resulting distribution function after the decoupling 
of our new light species, we then calculate the energy density at recombination, which is 
reported as the contribution to clS db function of the coupling parameters of the theory. 

As discussed in subsection III D\ each new light species must also decouple after the QCD 
phase transition in order to be constrained by Planck, which limits the dimensionality of the 
operators we choose to consider. If our new species couples to the SM with an operator of 
scaling dimension d, the operator is suppressed by A A ~ d , where A is the approximate cutoff 
scale of the EFT. Dimensional analysis then indicates that, given independent experimental 
constraints, only operators of dimension d < 6 will be able to maintain equilibrium between 
a new species and the SM until after the QCD phase transition. 

Finally, we discuss possible extensions to the SM which do not contain new light species, 
but instead alter the neutrino distribution, through such means as decay or neutrino asym- 
metry. These models then enhance the neutrino energy density relative to SM predictions, 
leading to an increase in g*. 

A. Spin-0: Goldstone Boson 

The first possibility for new light species is a spinless scalar boson. However, the mass 
of any new scalar particle is generically sensitive to quantum contributions resulting from 
interactions. While supersymmetry could potentially preserve the naturalness of scalar 
masses, the observed particle spectrum indicates that any couplings between the SM and 
new light scalars would mediate supersymmetry-breaking mass terms significant enough to 
require fine-tuning. The only viable symmetry which can protect the mass term of such light 
scalar bosons is then a shift symmetry, — > + e. This is precisely the symmetry present in 
the Goldstone modes of a spontaneously broken global symmetry. In the limit of an exact 
global symmetry, the mass of the corresponding Goldstone boson is restricted to be zero, 
with any quantum corrections forbidden by the symmetry. Even if the symmetry is inexact, 
the mass of the pseudo-Goldstone is proportional to the symmetry-breaking terms in the 
original Lagrangian, rather than the cutoff of the effective theory. We therefore restrict 



25 



ourselves to the study of Goldstone bosons, as other theories of light scalars are generically 
tuned and unnatural. 

One example of theories containing Goldstone bosons is the QCD axion model, which 
was originally motivated by the strong CP problem. The axion is then the Goldstone boson 
corresponding to a new Peccei-Quinn (PQ) global symmetry. While the original axion model 
with an electroweak PQ breaking scale has been observationally excluded, other variants on 
this model could result in viable contributions to g*. Axion-like species are also ubiquitous 
in string compactifications (the so-called 'String Axiverse') 61] . These particles arise as the 
imaginary parts of moduli of closed cycles in the compactification space. The imaginary 
parts then inherit a shift symmetry which arises from the gauge structure of the UV theory. 
Some fraction of these moduli will potentially be present in the low-energy theory, presenting 
another possibility for new light species. The real parts of string moduli do not concern us 
as they tend to be stabilized by fluxes in string compactifications, and when they are not, no 
symmetry protects them from receiving quantum corrections to their masses, pulling them 
out of our sub-eV range of interest. 

As we will discuss below, collider and astrophysical bounds on such particles make the 
simplest models too constrained to ever come into thermal equilibrium with the SM thermal 
bath. There are small corners in (flavor-dependent) parameter space which are still viable 
and in which they could in principle have a small but non-negligible impact on the effective 
number of relativistic degrees of freedom. We conclude that unless we are very lucky, the 
addition of a natural massless or near massless scalar will have, at best, a tiny impact on 
the CMB and thus would require significant advances in our ability to measure g*. 

The restriction to Goldstone bosons automatically reduces the possible couplings between 
a new scalar and the SM. There are two classes of interactions, both of which correspond 
to dimension-5 operators. The first corresponds to a derivative coupling between <ft and 
a SM fermion, such as the electron. This scalar could generically couple to only quarks, 
only leptons, or both quarks and leptons. When there are quark couplings in the UV, 
there will be induced couplings of the Goldstones to baryons and mesons after the QCD 
phase transition. We consider each of these scenarios. The second case corresponds to 
an interaction between <p and a SM gauge boson. In particular, since we only discuss 
thermodynamics after the QCD phase transition, the only available gauge boson coupling 
is to the photon. This loop-suppressed term is present only when the global symmetry is 
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TABLE I: Charges of fields in the Goldstone boson model. and <3? are new scalars in the theory. 
h is what gives rise to the Higgs field in the low energy theory, and tp acts as one generation of SM 
leptons. 

anomalous. In the low energy effective theory, therefore, any combination of lepton, pion and 
photon couplings may conceivably be allowed. We explore all of these possibilities, finding 
that current experimental constraints are such that almost all scenarios are excluded. The 
only remaining viable contribution of Goldstone bosons to g* arises in the specific model 
in which the coupling of to muons is much stronger than with other SM species and a 
model-dependent mechanism determines the <ft distribution at temperatures above the muon 
annihilation. 

Due to the shift symmetry, any coupling between an exact Goldstone boson and SM 
fermions must only contain derivatives of the field <fi. As a simple example of how such 
interactions arise, we consider the case of an additional global U(l) symmetry, which is then 
spontaneously broken. While there are many other possible theories with Goldstone bosons, 
the basic form of the resulting EFT is insensitive to those details and easily extended to 
the case of larger global symmetries. As such, the specific theory presented here is intended 
solely for illustrative and motivational purposes. The resulting EFT calculations are far 
more general, applying to a large class of Goldstone theories. 

In this simple UV completion, we add two new scalars $ and to the SM. The gauge 
and global charges of our particle content are shown in Table [H In this example, the SM 
field ip is one generation of lepton, such as electrons, but this model can be generalized in a 
straightforward manner to contain multiple generations and quarks. 

The addition of a new global U(l) now prevents the couplings of the SM fermions to the 
Higgs boson. The Lagrangian terms relevant to our discussion are then 
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- y&^L^ c R + ^ f 6$ + h.c. - M||0| 2 - V($) - V(/i). 

The mass Me corresponds to some scale above those of the SM, such that G can be 
integrated out of the theory to obtain 



If the potential V($) is such that ($) = ^, the global £7(1) symmetry will be spon- 
taneously broken. We can then parametrize perturbations about this symmetry-breaking 
ground state as $ ~ + p/y/2)e 2l< ^^ A , where p is the radial mode and (f) is the azimuthal 
mode. Written in terms of these new degrees of freedom, the Lagrangian is of the form 



uA 1 ( 25 ) 

- »^ /A ^ + ^- - 2 M > 2 - 

where the resulting mass term M p ~ A. Finally, we can make the field redefinitions ipi — >■ 

ammg 



£ => ^m 2 + ^^) 2 + i^i 2 + + ^iJp^R + 



, ./'A ■ 1 (26) 

2A 



Written in this form, we see that <fi acts precisely as the Goldstone mode of the broken U(l) 
and couples derivatively to the SM fermions. We have also induced the Yukawa couplings of 



the Higgs to SM fermions y' = y^r- Lastly, the additional massive field p will be removed 



from the EFT at energies far below M p . 



Using identities found in 



62|, we can also write our Lagrangian in Dirac notation, result- 
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ing in 



C d -(d^) 2 + {d^hf + - -^-*y 7 5 * - y7i f ** + /i.e. - (27) 

In this form, it is simple to see that the interaction is specifically a derivative coupling 
between and the axial current of One might suspect that some theories could potentially 
generate a similar coupling between and the vector current for However, any interaction 
of that form must vanish due to vector current conservation. The conservation of the axial 
current is broken by the mass term for meaning that the axial coupling does not similarly 
vanish. However, this does imply that any interaction rate involving the axial coupling is 
necessarily proportional to the fermion mass m, and thus vanishes in the m — >■ limit. 

We also notice that in this example model, the suppression scale A is set by the original 
symmetry-breaking expectation value, and is insensitive to any original Yukawa coupling 
of the theory. The couplings of to the SM would then be flavor-blind in this scenario. 
More complicated global symmetries or charge assignments could potentially result in flavor- 
specific couplings. However, a flavor-specific basis generically leads to interactions which mix 
generations. There are greatly restrictive constraints coming from flavor physics, as we shall 
discuss briefly below. 

In the low-energy effective field theory, there are three possible classes of couplings of 
the Goldstone to SM fields: couplings to leptons, pions, and photons, as mentioned earlier. 
Lepton couplings have the property that the interaction rate scales as ~ T well above 
the mass of the lepton. Consequently, T^/H increases with time instead of decreasing with 
time, and the interactions come into equilibrium, falling back out of equilibrium only as 
the leptons annihilate away and their number density drops. Pion couplings are potentially 
induced by couplings of the Goldstones to the quarks in the UV, and photon couplings 
are generically induced when the original global symmetry is anomalous. These rates scale 
as ~ T 3 , and therefore these interactions fall out of equilibrium at some temperature. 
Depending on the particle content and couplings of the UV theory, any combination of these 
couplings could conceivably be generated in the low-energy theory, and so we analyze all 
such possible cases below, summarizing our findings at the end of this subsection. 

We will first consider the interactions of <fi with SM leptons. Whether these couplings 
are present or not in the low-energy theory depend on whether the lepton couplings are 
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FIG. 4: Dominant interaction process for the Goldstone-lepton coupling. 

present in the UV. In the case where such UV couplings are present, we can easily evaluate 
the low-energy theory in a similar fashion to above, where the low-energy theory includes 
0, electrons, muons, neutrinos, and photons. Due to the A suppression of the derivative 
couplings, the interaction rate between and the SM will be dominated by processes which 
only involve one Goldstone interaction term, shown in Figure HI Note that, as this dominant 
process involves the emission/absorption of a photon, the interaction rate has no depen- 
dence on the coupling between and neutrinos. Because of this, the only relevant lepton 
interactions for are those with electrons and muons. 

In the relativistic limit, dimensional analysis would expect the interaction rate between 
and a SM lepton ^ to scale as ~ S. However, the broken axial symmetry for ^ 
restricts the interaction rate to take the form ~ ^p^; furthering the argument for ignoring 
neutrino couplings. Since is linear in T, the expansion rate will drop more quickly than 
the interaction rate as the universe expands. This implies that could have been out of 
thermal equilibrium after the time of global symmetry breaking, and come back into thermal 
equilibrium with the SM leptons at some point before the leptons annihilated. The details 
of lepton-0 interactions near the cutoff of the EFT depend on details of the UV theory, and 
so consequently we cannot make any predictions about the thermal interactions of and 
the SM leptons until recouples to the SM leptons (if at all) when ~ H . 

However, this treatment so far has only applied to the high temperature limit, with the 
assumption that T 3> m. As the universe cools to temperatures comparable to the relevant 
lepton mass, this simplified form for the interaction rate will be substantially modified and 
needs to be computed numerically. The interaction rate will begin to drop rapidly as the 
leptons annihilate away, redistributing their entropy amongst the remaining coupled species. 
If A is very large and substantially suppresses T^, will not have recoupled by the time 
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the leptons annihilate, meaning that <p will forever remain out of equilibrium. Additional 
couplings beyond the lepton-only couplings we consider in this paragraph would be needed 
in order to have SM-0 interactions. For each lepton, there is then some maximum A for 
recoupling. Any Goldstone boson with a larger A will not be reheated by the entropy 
redistribution and therefore cannot substantially contribute to g* at recombination. 

The distribution function for <p prior to recoupling is dependent on the original process 
of decoupling at high energies, which is then sensitive to details of the UV theory, including 
the relative timings of global symmetry breaking and inflation. It is then impossible to make 
fully model-independent predictions for the contribution of <p to g* in the case where <fr only 
couples to leptons. However, in a large class of models, including our simple motivational 
example, <fi will also couple to quarks and photons, which results in qualitatively different 
evolution. 

Since we consider temperatures below the QCD phase transition, we must examine the 
interactions between and mesons, specifically pions and charged kaons, resulting from 
quark couplings. At these temperatures, however, the number density of kaons will be 
much lower than that of pions, such that any 0-kaon interactions will be subdominant. We 



can then simply focus on those couplings which involve pions. Following [63|, the original 
coupling of cf) to quarks can be rewritten in terms of the axial quark current. After the 
phase transition to mesons, interactions with the quark current are replaced by those with 
the axial pion current, which is safe from QCD renormalization effects. The full Lagrangian 
can be expanded to leading order in A and f n and subsequently studied, and depends on the 
details of the flavor structure in the UV. We assume flavor-blind couplings, as flavor-specific 
couplings in the UV do not alter our predictions for the thermodynamic properties of </>, 
but such couplings must obey additional constraints coming from flavor physics. With this 
assumption, we find that those terms which dominate the interaction rate between <f> and 
pions are 

C D -^-n+n-d^rt + -^rOjr+^^V- + J^^-Q^^ . (28) 
oj n A oj^A oj^A 

We have defined the ratio r m = ™ d d +™ u , where m u and are the up and down current- 
quark masses. We use the approximate value r m = 1/3, based on lattice QCD calculations 



64j, as well as the convention f n = 93 MeV. Interactions of this type will potentially keep the 
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Goldstone boson in thermal equilibrium until the pions fully annihilate and redistribute 
their equilibrium, depending on the suppression scale A. The corresponding interaction 
rates decrease more rapidly than the expansion rate, leading to the freezing out of the 0-7r 
interactions. 

Finally, there can be couplings of to photons via operators of the form 

c " -sAx;***"*- < 29) 

This operator arises because the axial symmetry in question can be anomalous. This 
A 7 is not necessarily precisely the same as the A which couples to SM fermions, though 
their orders of magnitude are similar in a large number of UV completions. This is because 
the operator can be induced by loops of SM fermions. The additional loop factor in the 
parameterization of A 7 is present because in these cases, the operator appears in the La- 
grangian suppressed by a loop factor relative to the fermion couplings. As mentioned earlier, 
depending on the UV structure of the model, this operator may or may not be present in 
the low-energy theory. Similar to pion couplings, this operator gives rise to a rate such that 
0-7 interactions to freeze out as we go to lower temperatures. 

We now outline the constraints on these scenarios, working in a general framework with no 
assumptions regarding the operator or flavor structure of couplings in the UV. The bounds 
are best stated in terms of the effective operators 

"-^'-^^ (30) 

where tyf can either be a charged lepton or the proton. The strongest bounds for these 
models come primarily from observations of stellar and supernova cooling, which will also 
greatly constrain other models within this work. The production of new light species which 
interact weakly enough to escape the interior of a star provides an efficient energy loss 
mechanism, affecting both stellar cooling and evolution. Comparison of SM predictions to 
astrophysical observations then provides a strong constraint on the interactions of such new 
species. The resulting constraints for Goldstone interactions are 
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A e > 2.9 x 10 9 GeV, 

A p > 3.5 x 10 9 GeV, (31) 



A 7 > 1.2 x 10 7 GeV. 



More details about these bounds can be found in 6514681] . The relationship between the 
effective proton scale A p and the UV quark coupling scale A = A q present in equation f[2"gj) 
depends on phenomenological parameters in the baryon chiral Lagrangian, as well as details 
of the UV theory. However, A p and A q are related by an 0(1) number. Consequently, we 
use the conservative bound A q > 5 x 10 8 GeV. 

In addition, there are constraints on the set of off-diagonal operators schematically of the 
form -£—d fl <f)fi'y fl 'y 5 e coming from /i — > e + $ js^]. These bounds restrict 

fie 

A Me > 1.6 x 10 9 GeV. (32) 

These operators' contributions to early universe thermodynamics are not significantly 
different from that of muon couplings during the era following the QCD phase transition. 
Consequently, we do not consider this case to be qualitatively distinct from the case with 
muon couplings, but considerably more constrained, and so we do not consider these oper- 
ators further. 

Finally, there are no direct constraints on A p , but there is an indirect constraint coming 
from the <fi contribution to the anomalous magnetic moment of the muon, which is induced 
by the interactions we consider at one loop. We estimate that the contribution of an axion 
to a M is 

2 

Sa u . (33) 



Current measurements on the muon anomalous magnetic moment 64j restrict A M > 350 
GeV. 

In order to consider the general list of all possible models, we present our results for 
each interaction separately. For a large number of models, multiple such interactions will be 
present, such that these results will be even more restrictive. 
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Electrons: In the case where there are no couplings to pions or photons, the distribution 
function for (f) immediately after the QCD phase transition is set by details of the UV theory. 
In particular, it depends on the details of global symmetry breaking, and also potentially 
on the details of inflation, depending on the relative timing of these two events. However, 
under the mild assumption that the Goldstone boson energy density p is not substantially 
larger than p$M after the QCD phase transition, we find that any interactions between 
and the electron which are strong enough to partially thermalize <fi require values of 
A e considerably lower than the current experimental bound. Therefore, we conclude that 
electron-0 interactions did not play a role in early universe thermodynamics. 

Photons: We find that the bounds on photon- Goldstone couplings are restrictive enough 
such that the Goldstone bosons could not have been in thermal equilibrium with the SM 
after the QCD phase transition if couples only to photons. If there are additional interac- 
tions with other SM species, the combination of loop suppression factors with the restrictive 
bounds on Goldstone-pion couplings and Goldstone-electron couplings are such that interac- 
tions with photons are too weak to play a role in post-QCD phase transition thermalization 
of the Goldstones. Therefore, regardless of whether such photon interactions are present or 
not, they are irrelevant. 

Pions: The presence of sufficiently strong pion-0 couplings allows for <p to be in thermal 
equilibrium with the SM after the QCD phase transition. In the case where the couplings 
to the pions are the only couplings present in the low-energy theory, we find that in order 
for Goldstone bosons to be in thermal equilibrium with the pions and receive any of the 
pion's entropy, the coupling suppression scale must be A < 5 x 10 6 GeV. This is illustrated 
in Figure [5j The maximum possible A necessary is far below the bound on A quoted above, 
and therefore the decoupling of the Goldstone must have happened during or before the 
QCD phase transition, making the Goldstone not a viable candidate for a contribution to 
g* in theories containing only pion interactions. 

Muons: In the case of muon-only couplings, like the case of electron-only couplings, it 
is not possible to give well-defined initial conditions for the Goldstone boson distribution 
function just prior to the recoupling of the Goldstone boson to muons, preventing us from 
making any model-independent statements about the thermodynamics of the Goldstone 
bosons during recoupling. However, it is possible to state upper bounds for Ag*. For all 
reasonable initial configurations of the Goldstone distribution function, if A M < 10 6 GeV, 
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Effective Scale A (GeV) 

FIG. 5: Ag* due to a single Goldstone boson which interacts with only pions. The contribution to 
5* at recombination is given as a function of the effective scale A, which suppresses this interaction. 
The gray region for A > 5 x I0 6 GeV corresponds to models which decouple during the QCD phase 
transition. The provided values of Ag* should therefore only be interpreted qualitatively in that 
region. Supernova and star cooling constraints on this scenario limit A > I0 9 GeV, and so this plot 
demonstrates that the Goldstone must have decoupled during or before the QCD phase transition. 

then we expect that the Goldstone will fully thermalize with muons, and receive the full 
share of the entropy redistribution, leading to a contribution of Ag* = 0.26, or AN e ff = 0.57. 
We plan to pursue more precise predictions in future work. 

However, if these couplings are present in conjunction with couplings to at least pions, 
then it is possible to study the decoupling of Goldstones from the SM, as the Goldstones had 
been in thermal equilibrium in the era leading up to muon annihilation. Note that electron 
and photon interactions with are too suppressed to play a role in thermodynamics in the 
era immediately following the QCD phase transition, and so the below results apply to any 
case where at least pion and muon couplings are present. In order for a Goldstone to have 
received any entropy at all from SM annihilations following the QCD phase transition, it 
must have coupled with A < 1.5 x 10 7 GeV. This is illustrated in Figure |6j This is also 
below the bounds quoted above, and therefore the Goldstone is not a viable candidate for a 
contribution to g* in this scenario. 



35 



0.25 - 



0.225 - 



bo 
<l 




0.175 - 



0.15 - 



0.125 - 



5xl0 5 lxlO 6 2xl0 6 5xl0 6 

Effective Scale A (GeV) 



1x10' 



2x10' 



FIG. 6: Ag* due to a single Goldstone boson which interacts with at least pions and muons. The 
contribution to at recombination is given as a function of the effective scale A, which suppresses 
this interaction. The blue-gray region for A > 1.5 x 10 7 GeV corresponds to models which decouple 
during the QCD phase transition. The provided values of Ag* should therefore only be interpreted 
qualitatively in that region. Supernova and star cooling constraints on this scenario limit A > 10 9 
GeV, and so this plot demonstrates that the Goldstone must have decoupled during or before the 
QCD phase transition. 

In the particular case of a QCD axion with flavor-blind couplings and A ~ 10 9 GeV, the 
couplings of the axion were weak enough that the axion must have decoupled before the QCD 
phase transition, where the number of degrees of freedom was at least 61.75. Consequently, 
it could not have contributed more than Ag* « 0.03, corresponding to a AN e ff « 0.06. 

To summarize, there are no parts of the minimal, natural parameter space where the 
Goldstones had been in thermal equilibrium with the SM through the QCD phase transi- 
tion, leading to <fi decoupling from the SM after the QCD phase transition, which do not 
directly conflict with bounds coming from star and supernova cooling. As such, the effects of 
Goldstone bosons on the CMB in the predictive part of the parameter space are well below 
the sensitivity of the Planck satellite. One can, however, have couplings to only the muon 
with A < 10 6 GeV, and obtain a nontrivial contribution to Ag*. However, there is no pre- 
dictive power in such a scenario, as it is not possible to give well-defined, model-independent 
initial conditions for the Goldstone distribution function. Therefore, the only viable set of 
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theories must contain a highly specific hierarchy of couplings, such that interactions with 
muons are much stronger than those with other SM fields present after the QCD phase 
transition, without the generation of significant off-diagonal couplings. 

B. Spin-i; Light Fermion 

Natural models of light spin-| fermions are made more easily than those containing light 
scalar bosons. This naturalness can arise due to chiral symmetry, which corresponds to a 
rotation of the field by an arbitrary phase, x ~^ eia X- This symmetry permits any fermion 
gauge and kinetic terms, but forbids Majorana mass terms. Even if chiral symmetry is 
explicitly broken by the presence of a small fermion mass, corrections to this mass parameter 
are in general proportional to the original value, eliminating the need for any fine-tuning. 
Similarly, Dirac mass terms can be protected by an axial symmetry. 

Because of this protective symmetry, there are many allowed interactions for light 
fermions. The possible models include interactions with SM gauge bosons, either through 
direct gauge couplings or dipole moments, as well as interactions with SM fermions through 
effectively pointlike operators, which result from the exchange of heavy intermediary parti- 
cles. 

While the focus of this work is on model-independent predictions using EFT couplings 
of new light species, we also consider the possible contribution to g* of two well-motivated 
specific models of new light fermions. The first corresponds to the addition of two light 
sterile neutrinos to the SM, consistent with results from short baseline neutrino oscillation 
experiments. The second corresponds to the supersymmetric partner of the axion, the axino, 
which arises in a large number of supersymmetric theories. 

1. Gauge Interactions 

One possibility is that a new light fermion x is charged under the SM gauge groups. The 
coupling strength of a fermion in any representation of SU(3)c or SU(2) L is completely 
fixed by the representation theory of these groups. While x could naively have any value 
of hypercharge, the prospect of gauge unification indicates that hypercharge values are 
also discrete. Any new fermion in non-trivial representations of the SM gauge groups will 
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therefore couple with the same strength as the SM fermions. Light species which possess 
electromagnetic or color charge are completely excluded. The only remaining option is a 
neutral fermion, which must couple to the Z, but these light fermions are excluded by 



measurements of the Z- width 64]. As such, light fermions in any non-trivial representation 
of the SM gauge groups are excluded as potential candidates for contributions to g*. 

However, if x instead coupled to some new gauge boson, kinetic mixing between this new 
field and the SM gauge bosons would lead to mixing-suppressed SM gauge couplings for x- 
Any such 'millicharged' light fermion therefore requires the existence of a new gauge boson, 
which would also contribute to g*. We consider the details of new gauge bosons and the 
resulting millicharged interactions in subsection IIII CI 

2. Dipole Moments 

While a new fermion cannot carry SM charges, x could still interact via dimension-5 
dipole moment operators. The only nontrivial dipole interactions between x an d SM gauge 
bosons are those with the hypercharge gauge boson, which are of the form 

C^-jB^XL^XR + h.c, (34) 

where the structure of these operators is such that we must introduce two new Weyl fermions, 
Xl and Xr- These interactions can arise from loops involving heavy charged intermediaries, 
whose mass and couplings set the dipole moment scale A. 

However, the charged intermediary loops that generate this operator necessarily preserve 
only the vector £7(1) global symmetry of x-, which is precisely the symmetry structure al- 
lowed by a Dirac mass term mxLX C R- Therefore, any UV completion reducing to the theory 
containing the Lagrangian terms of equation (134]) must also allow for a Dirac mass term. It 
is not apparent how to create a UV completion of this model which induces only a dipole 
term corresponding to large mass scales, while generating the Dirac mass < eV in a natu- 



ral fashion. Experimental constraints from star cooling observations 69] currently limit x 
dipole moments to 

A > 10 9 GeV. (35) 



38 



Due to the resulting large separation of scales in this highly constrained EFT, we do not 
consider a theory with new light species possessing a SM dipole moment to be a viable, 
natural candidate for a contribution to Ag*. 



3. Four-Fermion Interactions 



Another possibility for EFT interactions of light fermions is the dimension-6 couplings 
of a single Weyl fermion x or a Dirac pair of fermions X to SM fermions. Such couplings 
can arise due to the exchange of a massive scalar or vector boson. Spontaneously broken 
gauge symmetries, which generate such massive interactions, are present in a large class of 
theories. 

As we will see below, light Weyl fermions with dimension-6 couplings are strong candi- 
dates for significant contributions to g*. Interactions suppressed by scales A ~ 2 TeV will 
keep new species in equilibrium until after the QCD phase transition, leaving such species 
with a detectable energy density at recombination. The strongest independent bounds on 
such models are placed by collider experiments, which will continue to probe the relevant 
parameter space. These theories will then potentially be discovered or fully excluded with 
the LHC. 

In Dirac notation, the possible dimension-6 operators present after electroweak symmetry 
breaking (EWSB) take four forms, 



^XX$$ (Scalar), 

X 

— X7 5 X^7 5 ^ (Pseudoscalar), 

f (36) 
— X^X^* (Vector), 

-^X 7 ^7 5 X^7 5 ^ (Axial), 

where \I/ corresponds to any SM fermion and the suppression scale A arises from the mass 
and couplings of the exchanged intermediary. We instead discuss the couplings of a Weyl 
fermion x below, as our results can simply be scaled by a factor of 2 to account for the two 
fermions in the Dirac case. As all four operators are dimension-6, the interaction rate will 
drop more quickly than H, leading to the decoupling of x from the SM as the universe cools. 
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FIG. 7: Ag* due to a single Weyl fermion which interacts with the SM via the exchange of a massive 
vector boson. The contribution to g* at recombination is given as a function of the effective scale 
A, which suppresses this interaction. Constraints on this operator are given for interactions with 
electrons (purple) and quarks (blue), which come from the LEP and LHC collision experiments. 
The green band indicates couplings which are 95% excluded by a Planck result of = 3.50 ±0.12. 
The gray region for A > 5 TeV corresponds to models which decouple during the QCD phase 
transition. The provided values of Ag* should therefore only be interpreted qualitatively in that 
region. The results for scalar, pseudoscalar and axial couplings are effectively the same. The results 
for a Dirac fermion are double those given in this figure, indicating that they must have decoupled 
during or before the QCD phase transition to be compatible with the Planck data. 

Couplings of this new fermion to quarks can induce couplings to pions in the low-energy 
theory. However, any interaction arising from the scalar or pseudoscalar operators will 
not be protected against strong renormalization effects, such that we cannot make precise 
theoretical predictions for Ag*. If such a species is independently discovered, potentially in 
collider experiments, and these interactions are precisely determined, a detailed calculation 
could then be performed. In addition, the vector or axial interactions are such that mesons 
will have no charge under such couplings, with no induced pion couplings in the EFT 7 . 

7 Vector or axial interactions between pions and \ would result from models with couplings which are not 
flavor-blind. Such couplings can only arise from the spontaneous breaking of nonabelian gauge groups 
which do not commute with flavor symmetry. Such models require a significantly larger particle content, 
thus violating our minimality requirement. 
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For each of these interactions, therefore, we can scan over possible effective suppression 
scales. The resulting contribution to g* as a function of A is given in Figure [7J for only the 
vector coupling, as the results for all four models are equivalent to within 5%. Therefore, any 
distinction between these models is below the experimental resolution of Planck. We also 
assume identical couplings to electrons and muons. For the possible case of flavor-specific 
couplings, the resulting Ag* will be the same as the flavor-blind case, where the equivalent 
flavor-blind A is the smallest flavor-specific A. 

The strongest experimental constraints on such couplings are indicated in Figure [7J for 
both electron and quark interactions. These bounds come primarily from $ T + mono- 
jet /monophoton searches at LEP and the LHC, again assuming universal coupling to quarks. 
The LHC bounds specifically came from 10 fb _1 of data, so we expect these experimental 
suits to improve in the near future. Details of these exclusion limits can be found in 
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re- 



72|. 



Couplings to muons are largely unconstrained in a flavor-specific model, but in the universal 
coupling case, constraints on any species would therefore limit the muon interactions. 

As we see, theories with effective suppression scales A > 5 TeV decouple prior to the 
muon entropy redistribution and are therefore predominantly affected by the QCD phase 
transition. As such, our results beyond those scales can only place an approximate upper 
bound on the possible contribution Ag*. However, there is a range of potential suppression 
scales below 5 TeV but above the current experimental bound which is compatible with the 
constraints coming from a Planck measurement of g* = 3.50 ± 0.12. A model with a light 
Weyl fermion with dimension-6 interactions with SM fermions is therefore a viable model for 
substantial contributions to g*. Our results indicate that a Dirac fermion contributes double 
what a Weyl fermion does at the same A, and that scalar, pseudoscalar, vector and axial 
vector operators give the same results to within 5%. Consequently, Dirac fermions must 
have decoupled before or during the QCD phase transition in order to be compatible with 
the data from Planck. Future results from the LHC will continue to probe these interaction 
scales, providing an independent means of discovery or exclusion of such models. 

4- Sterile Neutrinos 

One well-motivated model of new light fermions is that of sterile neutrinos. The observa- 
tion of neutrino oscillations necessitates an extension to the SM. A possible explanation of 
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these oscillations is the addition of some number of sterile neutrinos v c Rl which are SM gauge 
singlets. In most models, the additional states beyond those of the SM are taken to be quite 
massive, such that they are absent from any low-energy theory. However, it is possible for 
these states to be much lighter and therefore potentially contribute to g* at recombination. 

One possible extension to this framework is the inclusion of a new massive gauge boson 
Z', correspon ding to a spontaneously broken U(l), which couples to both the SM and sterile 



neutrinos 
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761 ] . Gauge anomaly cancellations would require the presence of three sterile 



neutrinos, but the exact mass spectrum is not explicitly fixed. The Z' would be integrated 
out of the low-energy theory, generating the dimension-6 operators discussed in the previous 
section, such that those results could be applied to this model. If such a Z' were discovered 
at the LHC, with a mass such that sterile neutrinos decoupled during the QCD phase 
transition, this model would then provide an observational probe of that phase transition. 
Detailed analysis of this particular model is beyond the scope of this paper, but will be 
pursued in future work. Here, we focus only on the addition of sterile neutrinos which are 
SM singlets, with no additional gauge groups. 

Multiple short baseline oscillation results suggest the existence of neutrino mass splittings 
distinct from those required to fit solar and atmospheric neutrino data, perhaps requiring 
additional light neutrino states (for details, see |77j and references therein). Motivated by 
multiple analyses of this short baseline data, we consider the addition of two sterile neutrinos 
with masses ~ eV to the SM. Possible alterations to this particular framework can be found 
in 78M8 01] . While previous work has explored the decoupling and contribution to g* of sterile 



neutrinos in a limited framework 
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-851]. we apply our full numerical algorithm to find the 
precise evolution of this particular model. 

As we shall see, sterile neutrinos with masses near the scale of recombination will decouple 
from the SM near the end of the muon entropy redistribution. The resulting energy density 
of such states will then be comparable to those of active SM neutrinos. However, as these 
species will become nonrelativistic during recombination, they do not affect the CMB in 
precisely the same manner as massless particles. The contribution of sterile neutrinos to 
the CMB can therefore not be accurately parametrized solely in terms of the single number 
g*, and a more detailed analysis of their effect on Silk damping and the early ISW effect is 
necessary to discover or exclude these species with CMB measurements. While this approach 
will be pursued in future work, here we present a generic means of determining the energy 
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density and pressure of sterile neutrinos as a function of the SM temperature, which can 
potentially be mapped to the details of CMB formation. 

Sterile neutrinos can appear in the Lagrangian in the terms y l J u^hU^. After EWSB, these 
interactions give rise to Dirac masses between the v c R and the active left-handed neutrinos 
u Ll as well as interactions with the physical Higgs suppressed by the very small couplings 
y u . The simplest framework for distribution functions is the description of propagation 
eigenstates, which have a well-defined relation between momentum and energy. We there- 
fore diagonalize the propagation basis and discuss the couplings of the resulting neutrino 
propagation eigenstates. 

In the framework where there are no new gauge bosons and only a Dirac neutrino mass, 
the presence of sterile neutrinos does not give rise to a substantial Ag*, due to the fact that 
the left-handed neutrino state is a Hamiltonian eigenstate. As an illustrative example, 
we consider the case of one generation, though the same result arises for more neutrinos 
with a general flavor structure. The neutrino mass terms are therefore simply 

C D —mu c R UL + h.c, (37) 

and the only interactions between u c R and the SM are suppressed by the small neutrino 
coupling to the Higgs. Consequently, we work in the free theory to study the effects of 
neutrino oscillations. The mass eigenstates are 

k±> = ^(K>±K», (38) 

where {v r \vl) — 0- Since the Hamiltonian is H = ^p 2 + m 2 , and only the state \u L ) will 
be produced from electroweak interactions, we can ask about the amplitude for it to be 
measured as a \u c R ) after some time difference t. We obtain 



We- iH >L) = ^=K\e- iHt (\u + ) + \u.)) 

y/T Rl V 1+7 1 V (39) 

= e- i( v / P 2 + m2 ) t (^|z/ L ) 

= 0. 
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Consequently, the only interactions which maintain equilibrium with v c R are those which 
are greatly suppressed by the neutrino Yukawa couplings. These states then decouple from 
the SM at very high temperatures, leading to a small Ag*. 

While this is true for the case where there is only a Dirac mass, there is also no gauge 
symmetry which prohibits the existence of additional Majorana mass terms for the sterile 
neutrinos. After EWSB, the most general mass terms allowed for a theory with multiple 
generations are therefore 

C D -m i3 u c Ri u Lj - ^MijU%^ + h.c, (40) 

where the Dirac mass matrix is simply determined by the neutrino Yukawa coupings such 
that rriij = vy^. The Majorana masses Mj could be forbidden by enforcing a total lepton 



number symmetry, but we will not require this symmetry. Following the approach of 86- 



, we also only consider renormalizable mass terms, rather than including higher order 
Majorana masses for the active neutrinos. 

Again, the propagation basis must first be diagonalized. The resulting propagation eigen- 
states will then consist of linear combinations of gauge eigenstates, with couplings dependent 
on the original mass terms. We study this in complete generality below, but as a conceptual 
example, we again consider the case of one active neutrino and one sterile neutrino. There 
are only two possible mass terms, a Dirac mass m and Majorana mass M, with the relevant 
Lagrangian terms of the form 

C D v\i$v L + v R \$v c R - g z v\$v L - g w e\ffiv L - mv L v c R - ^rV R v R + h.c. (41) 

The mass matrix can be diagonalized by a simple rotation of vl, v c r into the new states z/ 1; 
v-i. However, the resulting mass for v\ is tachyonic, which then requires the field redefinition 
v\ —> iv\. The resulting new states are 



M 

V4771 2 + M 2 



M 

V4T71 2 + M 2 ' 



-p 



u 2 



V2 



M 

^Am 2 + M 2 



Vl + \ / 1 + 



M 



VAm 2 + M 2 



R > 



'R 



(42) 
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The corresponding Lagrangian terms are 



C D u 



1 



V4m 2 + M 2 - M) v x v x - - [V^m 2 + M 2 + M) v 2 v 2 



4 



■ 9w L . 
+ z— pWl + 

-fK 



V4m 2 + M 2 



4m 2 + M 2 



V4m 2 + M 2 



M 2 



M 



M 




v\$V2 — V % $V\ ) + h.c. 



(43) 



As mentioned earlier, the gauge couplings of each eigenstate are dependent on the original 
mass matrix. Based on experimentally-indicated mass splittings, we expect the approximate 
mass scales m < 0.1 eV and M ~ 1 eV. This slight hierarchy of scales allows us to gain 
some general insight by considering the so-called 'seesaw' approximation M > m, as the 
Majorana term is the dominant mass contribution. In the literature, it is most common to 
apply this seesaw mechanism to models with the Majorana mass M beyond any mass scale 
of the SM, such that there are no additional light degrees of freedom. However, even for 
M ~ eV, this approximation is still sufficient, but now results in the presence of new light 
species. 

To leading order in ^, we then obtain 



In this limit, there are a few important things to notice. From equation (|42j) . we see that 
the light species v\ is almost entirely composed of v^. The gauge couplings for this pre- 
dominantly active neutrino are therefore largely unsuppressed, while its mass is suppressed 
relative to the original Dirac mass m. The heavy species v 2 is mainly composed of v c R . This 




(44) 
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predominantly sterile neutrino has suppressed gauge couplings and a mass of approximately 
M. The dominant interactions of v 2 are those which involve only one sterile species, as these 
are only suppressed by This includes interactions between u 2 , charged leptons, and W, 
as well as those between u 2 , Ui, and Z. Such couplings introduce flavor-changing neutral 
currents into the SM, but these are not significantly constrainted in the neutrino sector. 

2 

Any couplings which involve two v 2 states are instead suppressed by j^. 

The general results of the seesaw approximation are therefore a very light, mostly active 
species v\ and a heavier, mostly sterile species v 2 . In the specific case of interest, where the 
new sterile state is still a light degree of freedom, we easily obtain the approximate masses 
mi ^ < 10~ 2 eV and m 2 ~ M ~ 1 eV, consistent with the experimentally-indicated 
mass splittings. 

We can now generalize this toy model to include the three SM active neutrinos and some 
number Nr of sterile neutrinos. The full mass matrix can be diagonalized by a rotation 
of the original neutrino gauge eigenstates by a unitary matrix U. Such a rotation induces 
couplings of all six propagation eigenstates to the Z. However, the resulting couplings are 
not at all constrained by measurements of the Z-width 64J], provided there is sufficient 
separation between the neutrino masses and the Z mass, as is the case in our model. To see 
this, we define the original Z-coupling matrix to be gzG, where G = diag(l, 1, 1, 0, 0, 0). The 
number of neutrinos as measured by the Z-width is then simply counted by Tr(G^G) = 3. 
After applying the Takagi diagonalization procedure to the mass matrix, G — > G' = UGU\ 
with U unitary. The new measured number of neutrinos is then Tr(G^G') = Tr (G^G) = 3. 
While the Z-width is therefore unaffected by the presence of additional light states, the 
thermal properties of the early universe are, such that light sterile neutrinos can contribute 
significantly to g*. 

If there are Nr sterile neutrinos added to the three active neutrinos of the SM, there are 
then 3 + Nr mass eigenstates. For Nr < 3, a generic set of original mass terms will result 
in 2Nr non-zero masses, 2Nr independent rotation angles, and 3(Nr — 1) CP- violating 
phases. Neutrino oscillation experiments have currently measured two mass splittings and 
three rotation angles, requiring Nr > 2 for this sterile neutrino model. 

Motivated by the recent short baseline results, we consider the effect of two sterile neu- 
trinos to Specifically, we use the experimental best fit parameters from 86]. While this 



particular experimental fit contains a third sterile neutrino with mass ~ keV, which then 
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serves as a potential warm dark matter candidate, that species can also be completely de- 
coupled from the low-energy theory by further raising its mass, without affecting the other 
parameters of the theory. There then remain two possibilities for the low-energy theory to 
fit the data, which correspond to either the normal (NH) or inverted (IH) hierarchy. The 
masses of the two predominantly sterile neutrinos are 0.7, 0.9 eV for NH and 0.8, 1.2 eV for 
IH. The remaining best fit parameters can be found in [86], specifically the NH1 and IH2 
models. The details of this model, specifically the CP violating phases, are unimportant for 
our results, as CP violation only arises at loop-level in this model. The relevant parameters 
are the neutrino masses and mixing angles, which are consistent with multiple analyses of 



the short baseline data |77l. I87H93]. 



The behavior of our simple example continues to apply, with the gauge couplings of 
the additional mass eigenstates suppressed by factors of \/ ^ m " ~ 0.1. These suppressed 
couplings cause the additional neutrinos to decouple from the SM prior to the decoupling 
of the predominantly active neutrinos. This decoupling occurs right at the end of the 
redistribution of entropy from muons among the remaining species, such that the energy 
density of each mostly sterile neutrino is comparable to that of one active neutrino. However, 
the distribution functions for the mostly sterile neutrinos contain deviations from the pure 
Fermi-Dirac form at low values of p. 

If this were the end of the story, we would then predict the contribution to g* to be such 
that AN e ff ~ 2. However, this particular model is more subtle than any model analyzed 
previously. Up to this point, we have considered any new species to be precisely massless. 
However, this is no longer true for the specific case of sterile neutrinos. As such, we have to 
consider the effects of nonzero masses on the measurement of g*, specifically the differences 
between Silk damping and early ISW. Both of these processes are sensitive to the precise 
evolution of H, whose time-dependence is altered differently by additional massive particles 
than by new massless particles. 

As an illustrative example, suppose we had a sterile neutrino which instantaneously de- 
coupled at a temperature T<j 3> m v , when the universe had a scale factor a^. The distribution 
function would have been purely Fermi-Dirac at this time, and would not have been sensitive 
to m v . Consequently, at a fixed temperature T<j, it takes the form 
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/(*,?)= (expQQ+l) = (exp(i-v^T^) +l) . (45) 

However, as discussed earlier, the Boltzmann equation for a species which is noninteract- 
ing demands that / must only be a function of q = ap. Therefore, the unique distribution 
function which satisfies the above boundary condition and solves the Boltzmann equation is 



f(t,p) = 




at all later scale factors a. From this distribution function, we can compute the energy 
density and pressure, in order to discuss the contribution of this species to the evolution 
of the universe. This is the strategy employed below, using instead the output distribution 
function of our code to compute p and P. 

Silk damping is primarily sensitive to the overall expansion rate, and therefore the overall 
energy density, near the point of recombination. These additional neutrinos will add more 
energy density than is predicted solely by the SM. However, due to their mass, the new 
neutrinos behave as relatively hot dark matter, as they will have become nonrelativistic by 
the modern era. Consequently, they contribute to measurements of QpMh 2 today, whereas 
they did not impact the CMB in the same fashion that that dark matter did. The values of 
fl VR h 2 today, as determined by the final distribution function in our code, are £l UR h 2 = 0.017 
for the NH, and Q VR h 2 = 0.021 for IH. Consequently, the sterile neutrinos cannot act as 
all of the dark matter today, but rather contribute a small but non- negligible fraction of it, 
indicating a need for an additional species which constitutes the majority of dark matter. 
The exact contribution of a massive species to Silk damping is therefore sensitive to the 
amount of dark matter in our universe, which is dominated by uncertainty in the overall 
dark matter content, and a more careful analysis of the effects of sterile neutrinos is needed. 

Similarly, the early ISW effect is sensitive to the radiation/matter ratio following the 
formation of the CMB. Massive sterile neutrinos will be transitioning to a nonrelativistic 
distribution during this period, behaving as neither pure radiation nor pure matter. The 



48 



contribution of these neutrinos to early ISW is determined by both their energy density and 
pressure, which are determined by their process of decoupling. Again, though, the exact 
prediction of early ISW effects is also dependent on the precise energy density of non-sterile 
neutrino dark matter. 

The main complication to the calculation of Ag* for sterile neutrinos arises from the use 
of the specific ACDM framework in calculating cosmological parameters from CMB data, in 
which the mass of dark matter is significantly higher than the temperature of recombination. 
This leads to model-dependence in the reported bounds, which do not necessarily exclude 
models which fall outside of this framework. It is important to stress that the difficulty 
arises due to uncertainty in the precise expansion rate and dark matter content, not due to 
any calculational uncertainty in the sterile neutrino sector. In order to fully probe the effect 
of sterile neutrinos on the CMB, a more general analysis of the CMB anisotropy data must 
be taken, which includes the possibility of nonzero masses for these additional light species. 
Such an analysis is beyond the scope of this paper, but will be pursued in future work. 

However, what can be done currently is a calculation of the energy density p and pressure 
P of these additional species as a function of the photon temperature T, which does not rely 
on the structure of the ACDM framework. This is determined by using the best-fit param- 
eters taken from jsf]] to calculate the sterile neutrino distributions following decoupling, as 
stated above, then evolving forward to arbitrary temperatures. These two functions com- 
pletely characterize the contribution of additional neutrinos to the evolution of H. Figure 
IE] provides both the actual values for p and 3P near recombination for these massive sterile 
neutrinos, as well as the values obtained if the species were massless, but the couplings to 
the SM and the distribution functions near decoupling were held fixed. 

5. Axinos 

In supersymmetric axion models, the axion is part of a larger supermultiplet, which also 
contains a Weyl fermion, the axino, and a real scalar, the saxion. If supersymmetry is pre- 
served, the Goldstone mass protection of the axion extends to the axino and saxion, keeping 
them naturally light. However, supersymmetry-breaking effects then spoil this protection, 
generically lifting the saxion out of the sub-eV range of interest. As there are many models 
in which the axino can remain light due to a chiral symmetry, we consider the potential for 
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FIG. 8: Energy density p and three times the pressure 3P of two sterile neutrinos as a function 
of the photon temperature, for the (a) normal and (b) inverted hierarchies. The effect of nonzero 
masses is shown by comparing the actual p (blue) and 3P (red) of massive neutrinos to those 
obtained if the sterile neutrinos we considered were massless (gray). At high temperatures, p and 
P approach the massless limit, such that p ~ 3P. At temperatures below the neutrino masses, p 
begins to fall less rapidly than the massless case, while P falls more rapidly. 

them to serve as a new light degree of freedom. 

Below the global symmetry-breaking scale A, the couplings of the axino a to the SM are 
determined by both the axion shift symmetry and supersymmetry. We first consider the 
lepton terms. The superpotential terms giving rise to viable couplings are 



W ~ e A / A y?E<rH d L 3 + e^Vtf^, (47) 
where A is the axion superfield containing a, and the remaining fields are SM superfields. 
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The second term can only arise in theories with R-parity violation. At low energies, following 
EWSB and the removal of all SM superpartners, we then obtain the interactions of the form 



C ~ - J^U^ty - ^~av L . (48) 
A 2 m~ A 

The first term corresponds to a dimension-6 interaction between the axino and leptons, 
generated by the exchange of a slepton with mass mj. However, since A is the same scale 
suppressing axion couplings to leptons, it is constrained to be A > 10 9 GeV. If axinos only 
interact with the SM via this first operator, they will decouple long before the QCD phase 
transition and therefore not contribute to g*. The case of couplings to quarks is constrained 
with similar bounds, with the same results. 

The second term corresponds to a Dirac mass between the axino and SM neutrinos. In 
this model, the axino therefore acts as a light sterile neutrino, with a potential contribution 
Ag*, as previously discussed. 



C. Spin-1: Gauge Boson 

A massless spin-1 particle has fewer degrees of freedom than a massive one, and thus 
perturbative quantum effects cannot generate a mass, rendering a massless spin-1 particle 
technically natural. These gauge bosons are then automatic candidates for new light species. 
While gauge bosons can potentially acquire a mass through the Higgs mechanism, masses at 
scales < eV are generically unnatural, unless there is a more complicated particle content 8 . 
However, such non-minimal solutions are beyond the scope of this work, so we assume that 
any additional vector bosons are precisely massless. Similar to the case of a Goldstone 
boson, the corresponding gauge structure automatically restricts the available interactions 
for light spin-1 particles. The only possible operators are direct gauge couplings or dipole 
moment interactions with SM fermions, as well as kinetic mixing with SM gauge bosons. 

8 It is, of course, possible to Higgs the group at the TeV-scale, but have such a small gauge coupling that 
its mass is sub-eV (g < 10~ 12 ). However, such a small gauge coupling implies that it will only recouple 
at very low temperatures, and even then, only to neutrinos. As neutrinos would have already decoupled 
from the SM, such interactions can only redistribute the neutrino energy density and cannot increase the 
total energy density. Thus there are no contributions of such a model to g* . 
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1. Kinetic Mixing and Gauge Interactions 



As we will show, new massless gauge bosons with renormalizable couplings to SM fermions 
are viable candidates for contributions to g*. However, such models must include both kinetic 
mixing between the new and SM gauge fields, as well as the presence of new fermions charged 
under the new gauge group. The additional fermions obtain millicharged couplings to SM 
gauge fields, with astrophysical constraints such that these fermions must have masses > 
MeV. Such models are sensitive to the details of the full UV theory, as the hidden sector must 
come into equilibrium with the SM after originally being completely decoupled. The class 
of viable models is then constrained to a particular region of model-dependent parameter 
space. 

For minimality, we consider the addition of a single new U(l) gauge boson A', with 
associated field strength A'^ u . Long-range force constraints greatly restrict the possible 
direct gauge couplings of SM fermions to A', such that this model cannot contribute to g* 



941 ] . However, it is possible for A' to kinetically mix with the hypercharge gauge boson B 



with the following operator 

C D —A'^B^ (49) 
where e is simply a dimensionless mixing parameter. Such hidde n sec tor U(l) gauge bosons 



95Ml01] 



which mix with hypercharge arise naturally in many models 

This term indicates that our originally defined fields A' and B are not propagation eigen- 
states, and must be redefined to diagonalize the propagation basis. If both gauge bosons are 
precisely massless, then there is always a linear combination of gauge fields which does not 
couple to the SM. We can always define this linear combination as A' and the orthogonal 
combination as B, such that A' does not couple to the SM. 

However, if A' originally interacts with some new fermion Xi an y field redefinition will 
generically result in couplings between \ an d the SM gauge bosons. The new fermion can 
then act as an intermediary between A' and the SM, keeping all species in equilibrium. 
The most minimal theory involving new direct gauge couplings must contain both a new 
gauge boson A' and a new fermion \- More complicated hidden sectors are possible, but we 
restrict ourselves to this simple model. The particle content and charge assignments of this 
example model are shown in Table HI1 To allow for the possibility of a mass term for x, we 
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TABLE II: Assigned charges for A' model 

have included two new Weyl fermions, xl and Xr- F° r illustrative purposes, we will focus 
on the couplings of the right-handed SM leptons ip R and the right-handed Xr- However, the 
remaining fermion couplings will follow the same behavior and can be easily calculated. 

Following EWSB, the mixing term in equation ( 149|) will lead to kinetic mixing of A' with 
the photon and Z . The relevant Lagrangian terms are then 



C D _ -F^F^ - -Z^Z^ - jA'^A'^ - "-cwA'^F^ + e -s w A'^Z„ v 

\ 4 4 2 1 (50) 

- -miZ„Z" + r R \i$ + e4- e-^$W R + xt{i$ + 9aA')Xr, 

where sw and cw are sine and cosine of the weak mixing angle, and is the new A' gauge 
coupling. A refers to the photon field, and F is the electromagnetic field strength. As we 
can see, the gauge bosons are currently not described by propagation eigenstates. However, 
the kinetic and mass terms can be simultaneously diagonalized via repeated rotations and 
rescalings of gauge fields, after which one can perform an 5*0(2) rotation between A and 
A'. Such rotations preserve the kinetic terms but alter the couplings between specific gauge 
bosons and specific fermions. Again, it is possible to choose a basis such that A' does not 
couple to the SM. To first order in e in this basis, we then find 



^- ~F^F, V - -Z^Z, V - jA'^A' - \m\Z,Z» 

sw , (51) 

c w 

We see that the SM fermions retain their original SM gauge couplings and do not interact 
with A', while x now interacts with both the photon and Z. Specifically, the interactions of 
X to the SM depend on both its coupling to A' and the original mixing between A' and B. 
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A shift in the ratio due to diagonalization only appears at (9(e ). 

This particular choice of basis is technically arbitrary. It is also possible to redefine gauge 
fields such that the SM fermions possess millicharged couplings to A' and x possesses no 
couplings to the photon. The physics must be and is independent of the choice of basis. 
These rotations do not affect any physical observable, provided the observable is phrased in 
a basis-independent manner. Thermodynamic observables such as the overall energy density 
of massless gauge bosons, and therefore their contribution to g*, are also basis- independent. 
We specifically choose to work in the basis of equation fl5Tl) . where A' does not interact 
with the SM, because the resulting early universe thermodynamics are more transparent. 
However, it is important to stress that the same results are true, but less obvious, in other 
bases. 

The dominant interactions between x an d the SM are dimension-4 gauge couplings with 
the photon, as any interactions with the Z are suppressed at temperatures below the weak 
scale. Dimensional analysis then implies that at temperatures large compared to m x , the 
interaction rate is linear in temperature, T x ~ T. Similar to the Goldstone couplings to 
leptons, this means that at high temperatures x wn l be fully decoupled from the SM and 
then potentially recouples as the universe cools and the expansion rate drops when e < 1. 
Unlike the Goldstone case, x is always interacting with A', provided the A' coupling is 
sufficiently large, such that x an d A' can maintain equilibrium distributions. Therefore, 
the hidden sector has a well-defined temperature. The precise ratio of temperature of the 
hidden sector to the temperature of the SM prior to the recoupling of the two sectors is 
model-dependent, as more complicated hidden sectors will generally result in a wide range 
of possible temperatures. Consequently, we choose to explore a wide range of such initial 
ratios. 

The thermodynamics are sensitive to whether A' and x are i n equilibrium, rather than 
the precise coupling qa, so we can simply fix the value of qa to be sufficiently large, without 
loss of generality. We select the value g\ = 0.1, but our final results can be simply related 
to other values of qa- Once qa is fixed, there are only three remaining parameters that can 
change: the kinetic mixing e, the new fermion mass m x , and the ratio of initial temperatures 
Thid/T. 

Multiple star and supernova cooling observations, as well as va rious collider results, place 
significant constraints on millicharged fermions (for details see 98l . ll02j |). Specifically, models 
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with m x < 100 keV are restricted to e < 10~ 13 , such that these species will never thermally 
couple to the SM for all reasonable initial values of T^d- Light millicharged fermions can 
therefore not directly contribute to g*, but more massive fermions can instead indirectly 
alter the CMB by maintaining equilibrium between the SM and A', which then contributes 
a nonnegligible Ag*. However, for this to occur, we need m x < 150 MeV, such that \ is still 
present below the QCD phase transition. This therefore limits us to a very narrow range 
of allowed masses m x for models of millicharged species which affect the CMB. For models 
of this type, \ must couple to the SM prior to or during its annihilation, otherwise the 
hidden sector will again never couple to the SM. This limits the possible values for e and 
Thid for any given mass m x , in addition to constraints placed by independent observational 
and experimental bounds. 

To illustrate the general behavior of these models, we consider four possible fermion 
masses within the allowed mass range. The corresponding results are shown in Figure [91 In 
each case, the millicharged fermion has mass m x > 10 MeV, which are unconstrained by 
star and supernova cooling observations and therefore have the largest available ranges for 
e and Thid- The lowest mass shown in Figure [9] is actually m x = 50 MeV, as the results 
are equivalent for masses between 10-50 MeV. For each of these cases, we scan over possible 
values for the mixing parameter e, as well as possible values for the original hidden sector 
temperature T hid when the SM temperature T = 200 MeV. We specifically consider T hid 
below the SM temperature T, assuming a minimal hidden sector model containing less 
particle content than the SM. The hidden sector will thus be colder due to fewer entropy 
redistributions. This procedure involved a modified version of the original code, the details 
of which can be found in Appendix Q3] 

As we see in Figure [9j there is a basic pattern to the dependence of Ag* on both e and 
T^d. For very small values of e, the hidden sector is never coupled to the SM, and A' receives 
all of the x entropy redistribution. The contribution to g* is then dependent solely on the 
energy available in the hidden sector. The energy density increases as the initial temperature 
increases relative to the SM temperature. Initial temperatures of Thid > jT are excluded by 
a Planck result of g* = 3.50 ± 0.12. 

For increasing values of e, the hidden sector begins to couple with the SM, until at large 
values the two sectors quickly become fully coupled, regardless of the initial temperature 
Thid- In this regime, the contribution of A' to g* is precisely that of a new gauge boson which 
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FIG. 9: Ag* due to a single gauge boson which couples to a new fermion with mixing-induced 
SM gauge couplings. The contribution to </* at recombination is given as a function of both the 
mixing parameter e and the hidden sector temperature Thid when the SM temperature T = 200 
MeV. Results are presented for (a) m x = 50 MeV, (b) m x = 75 MeV, (c) m x = 100 MeV, and 
(d) m x = 125 MeV. Blue and purple regions to the left of the black line are allowed by a Planck 
measurement of g* = 3.50 ± 0.12, although these regions were never in thermal equilibrium with 
the SM. Regions to the right of the black line are excluded by this result from Planck. 

is originally coupled to the SM then decouples before the electron annihilation, Ag* f» 0.5. 

In the transitional region from completely decoupled to completely coupled, the contri- 
bution rapidly climbs to Ag* w 1, then rapidly decreases to the fully coupled limit for large 
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e. This enhanced contribution to corresponds to a fortuitous combination of coupling 
and mass values, in which A' participates in the muon entropy redistribution but is able 
to receive all of the x entropy. This occurs because x briefly couples to the SM, sharing 
the muon entropy, then quickly decouples as the muon and x number densities begin to 
plummet, such that the SM receives none of the x entropy. The result is a superheated 
population of A' bosons, which contain a large fraction of the total energy density. 

While the majority of this behavior has been largely m x -independent, we do observe a 
slight decrease in the transitional region values of Ag* as the x mass increases. Larger 
fermion masses result in the hidden sector decoupling earlier from the SM, and therefore 
receiving less of the muon entropy. Finally, there are no major distinctions between m x ~ 20 
MeV and m x ~ 50 MeV, as these masses are proximate to neither the muon nor the electron 
mass. 

For initial hidden sector temperatures below ^, the behavior will be largely unchanged 
from the low-temperature results presented here. Theories with small mixing parameters 
will remain fully decoupled and contribute negligibly to g*, while theories with larger e values 
will rapidly reach equilibrium with the SM, such that their contribution Ag* is insensitive 
to the initial temperature. 

We find that for any value of the initial temperature, e is restricted to be < 10~ 8 when 
there is a Dirac fermion x with masses between 10 — 150 MeV, forcing the SM and the hidden 
sector to never have been in thermal equilibrium. In addition, a scenario with m x < 10 MeV 
is inconsistent with constraints from star and supernovae cooling, and one where m x > 150 
MeV causes the A's to decouple before or during the QCD phase transition. The result in 
the absence of x is the same result as obtained by raising the x mass and integrating it out 
of the theory; there exists a basis in which there are no couplings between the hidden sector 
and the SM, thus preventing the thermalization of A'. The scenario of new gauge bosons 
which mix with SM hypercharge is therefore further constrained by results from Planck. 

2. Dipole Moments 

While it is always possible to eliminate any mixing-induced renormalizable couplings 
between SM fermions and a new unbroken gauge boson A'^, there could generically still be 
higher-order nonrenormalizable couplings after integrating out A'-SM interaction mediators 
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TABLE III: Assigned charges for A' dipole model, a is an unconstrained real charge, X an d <p 
are BSM heavy Weyl fermion and complex scalar fields, respectively, which generate the dipole 
couplings of A' to the fermions. %' is a vector-like partner for X to enable it to have a large mass. 

in the full theory. These interactions will dominate A'-SM interactions in the absence of 
other new light charged particles. 

Astrophysical bounds on such dipole interactions make the simplest models too con- 
strained to decouple at sufficiently low temperatures. Similar to the case of Goldstone 
bosons, the only viable models which contribute to g* then contain highly flavor-dependent 
parameters. 

If the low-energy effective theory contains no light species charged under U(1)a', then 
the dominant interactions between A' and the SM are of the form 

where A ltxv is again the associated field-strength tensor and M is the mass scale associated 
with the heavy species integrated out of the theory. After EWSB, the expansion of the Higgs 
field about its expectation value v will lead to dipole moment interactions of the form 

CD- -^A'^r^l + h.c. - jA'^j^^L + h.c, (53) 

where we have now denned an effective dipole scale A = 

An example of a UV theory giving rise to the low-energy dipole terms described above is 



£ D -yl^XL ~ y x X R hXL - vW^Xr ~ M( X ' l Xl + XrXr) + h.c. - M 2 |0| 2 , (54) 
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FIG. 10: One of the several diagrams which gives rise to a dipole coupling of SM fermions to an A' 

where xl and x c r are new fermions and is a new scalar, all with masses M 3> v. As shown 
in Table IHIl these new species are charged under both the electroweak symmetry groups, as 
well as the new U{1)a'- The fermions x'l an d Xr do not couple to the SM except through 
gauge interactions, and are present solely to give a large mass to xl and Xr an d cancel 
anomalies. The remaining fermions ipi and tp R are the three generations of SM leptons, 
though this model can easily be extended to also include quarks. Interactions of the form 
of equation fl52|) can be generated at temperatures far below M upon integrating out the 
0, X-, an d x' fields, through diagrams such as the one in Figure [10J After EWSB, the SM 
fermions ip obtain a nondiagonal mass matrix, which needs to be diagonalized by different 
unitary transformations on the left and right handed fermions, L and R. Therefore, the final 
induced operator, F^ip^a^ip^, has a y iJ flavor structure which is not simply y{y J 3 , but 
rather (L ■ yx) % (ffl ■ y^ ■ Therefore, a full parameter space for flavor couplings is achievable 
in the low-energy theory, and theories with hierarchies between electron- A' couplings and 
muon-A' couplings are easily achievable. 

The induced dipole couplings for pions must involve a composite pion operator which 
is antisymmetric in its two Lorentz indices. If the quark dipole moments are flavor-blind, 
such that the up and down quarks have the same couplings to A', then all such antisym- 
metric operators vanish. If, instead, the dipole couplings are not flavor-blind, interactions 
between A' and pions will potentially appear. However, there is no symmetry protecting 
against renormalization of such operators. We then expect these pion interactions to be 
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strongly renormalized, thereby preventing us from making robust predictions about such 
contributions to g*. We therefore focus solely on the dipole couplings of A' to leptons. 
With the interaction Lagrangian of the form 

C D -j-* f a^ f A'^ (55) 
where \I/ / can either be an elementary lepton or a composite nucleon, we obtain the bounds 



A P > 2.0 x 10 10 GeV, 

(56) 



A P ,n > 9.8 x 10 9 GeV. 



The se b o unds again come from star and supernova cooling, and details can be found 



in 



65 



103 



104] . There are also constraints on off-diagonal couplings A^ e coming from 



H — » e + $ [103| , which limit such couplings to 



A Me > 2.3 x 10 9 GeV. (57) 

There are no direct constraints on fx- A' interactions, but we may again estimate that A' 
contributes 

2 

Tfl 

**> ~ ieM (58) 

to the anomalous magnetic moment of the muon. This leads to the approximate constraint 
A M > 350 GeV. 

We find that the electron-only and electron-muon off-diagonal coupling scenarios are 
constrained to decouple before or during the QCD phase transition, preventing A' from 
contributing to in this part of parameter space. However, when the coupling of A' to 
those species in the SM present after the QCD phase transition is dominated by its coupling 
to the muon, there is a relatively large allowed range for A M . Our results for muon-dominated 
couplings are shown in Figure [TTJ We find that A M < 10 7 GeV in order for the A' to remain 
coupled after the QCD phase transition, and consequently contribute to g*. This requires a 
significant hierarchy between the electron- A 1 coupling and the muon- A 1 coupling, but such 
a hierarchy is compatible with an MFV-like framework, as the hierarchy does not need to 
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FIG. 11: Ag* due to an A' which interacts primarily with muons. The contribution to g* at 
recombination is given as a function of the effective scale A, which suppresses this interaction. 
The gray region for A > 10 7 GeV corresponds to models which decouple during the QCD phase 
transition. The provided values of Ag* should therefore only be interpreted qualitatively in that 
region. There are no direct constraints on this scenario, and so new exclusion limits can be placed 
through the results of the Planck satellite. The green region corresponds to values of A excluded 
by a Planck result of g* = 3.50 ± 0.12. 

be much larger than y e /y^. Values of A M < 10 6 GeV are inconsistent with a Planck result 
of = 3.50 ± 0.12, providing novel constraints on such a scenario. 

D. Spin-|: Gravitino 

Any model of supergravity contains the gravitino, which is the unique elementary spin-| 
particle. If supersymmetry were unbroken, the gravitino would be precisely massless. In a 
method similar to that of gauge symmetries, the spontaneous breaking of supersymmetry 
gives rise to a massless fermion, the Goldstino, which then becomes the longitudinal mode 
of the gravitino. As a result, the gravitino acquires a mass m 3 / 2 ~ ^f-j, where F is generally 
the largest supersymmetry breaking scale squared in the theory. Similar to the axino, the 
gravitino can then remain a light degree of freedom for sufficiently low supersymmetry- 
breaking scales. 
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Naively, the gravitino would interact solely with gravitational strength and would there- 
fore decouple at very high temperatures. However, at energy scales far above the gravitino 
mass, the Goldstino equivalence theorem ensures that the longitudinal components of the 
gravitino interact with Goldstino strength, pote ntia lly maintaining equilibrium with the SM 



105| , the Goldstino couplings to the SM are 



down to lower temperatures. As is well known 
of the form 

C D - JLxV^x^ (59) 

where T^ v is the stress-energy tensor comprised of SM fields. While this coupling is no 
longer gravitationally suppressed, it is still a dimension-8 operator, such that the gravitino 
will still decouple above the QCD phase transition for all viable supersymmetry-breaking 
parameters F and not contribute significantly to g*. 



E. Spin-2: Graviton 

The unique elementary spin-2 particle is the graviton. The graviton interacts solely with 
gravitational strength, such that it either decouples from the SM at very high temperatures 
or is never even in thermal equilibrium. Similar to the discussion of subsection III D\ gravitons 
which decouple at such large temperatures cannot contribute more than around 0.02 to 
well below the sensitivity of Planck. 



F. Models Without New Light Species 

Up to this point, we have considered the addition of new light species to the SM in 
order to increase the total relativistic energy density, p re i, at recombination. The other 
possibility is a modification of the distribution functions of light species already present 
in the SM. The distribution function of photons is well-established as Bose-Einstein by 
measurements of the CMB, such that the energy density of photons at recombination is 
known to high precision. The only remaining option is therefore a modification of the 
distribution function of neutrinos. If this change were to occur while neutrinos were still in 
equilibrium with the SM, then interactions with SM species would thermalize the distribution 
functions, washing out any original alteration. Therefore, new physics must only affect the 
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neutrino energy density at temperatures below the MeV scale. Here we briefly discuss the 
possible mechanisms which can alter the distribution function of neutrinos to increase 
at recombination: a neutrino asymmetry, decays of a new species, and new interactions 
between neutrinos and the remaining SM species. 

A simple modification to the distribution function of a species is the introduction of a 
chemical potential fi, such that 

/fe g )= e(E -, ) /T + 1 . (6°) 

for fermions, with fi — > —\i for antifermions. Such a chemical potential results in an asym- 
metry between the number of particles and antiparticles. Once a species has fully decoupled 
and freely evolves, the Boltzmann equation constrains / to remain solely a function of a(t)p, 
such that £ = ^ is then time- independent. We can then express any resulting effects in 
terms of this constant £. 

Although the neutrino distribution function is no longer Fermi-Dirac after decoupling 
from the SM, we shall assume it is for illustrative purposes. The total energy density stored 
in neutrinos and antineutrinos with nonzero £ is given by 



P» = (W-e 5 ) + Li 4 (-e «)) = 12Q ^1 + — + 



15? 



(61) 



7vr 4 

The presence of a nontrivial chemical potential for neutrinos would therefore increase the 
energy density, thereby increasing g*. The electron neutrino chemical potential affects the 
neutron-to-proton ratio prior to the start of BBN through reactions of the form p+u — > n+e. 
This ratio directly affects the helium-4 abundance after BBN, and so bounds can be placed on 
the electron neutrino chemical potential. Furthermore, since all neutrino mass eigenstates 
contain some wavefunction overlap with the electron neutrino, all of the neutrino mass 
eigenstate chemic al po tentials are constrained. The result is £ ^$ 0.1 for each of the three 



neutrino species 



106| and references therein). 



A second possibility is the late decay of some massive species to neutrinos [1071 ] . This 
heavy species must have fully decoupled at some higher temperature, leaving a significant 
relic energy density, with a decay rate such that it decays predominantly to neutrinos after 
neutrino decoupling but prior to recombination. These decays would then significantly alter 
the neutrino distribution, creating a large number of neutrinos with energies comparable to 



63 



the particle mass. The mass, number density, and coupling to neutrinos for this new species 
determine the precise contribution to g*, making this scenario highly model-dependent. 

Finally, the existence of higher-dimensional operators coupling the neutrinos to other SM 
species could potentially maintain thermal equilibrium between neutrinos and the SM until 
lower temperatures. A later point of neutrino decoupling would result in a larger share of the 
electron entropy being distributed to neutrinos, raising their energy density. The possible 
interactions with the lowest dimensionality are electromagnetic dipole moments or four- 
fermion interactions between neutrinos and electrons, which are significantly constrained by 
star cooling [69] and the LEP collider 72]. 



IV. CONCLUSIONS 



The Standard Model of particle physics represents our current knowledge of the quantum 
field theory that best describes all short-distance interactions down to 10~ 17 cm. Know- 
ing that this model is incomplete leads us to search for fundamental particles outside the 
Standard Model. While the search for heavier particles continues at colliders, we focus on 
another class of new physics - light, stable particles - which can be probed via their effects 
on cosmology, most strikingly on the Cosmic Microwave Background. In this article, we have 
surveyed what we call the most 'natural' (or least contrived) models and their parameter 
spaces. By doing so we lay out the reach of current and future experiments detailing the 
power spectra in the Cosmic Microwave Background and other probes of the initial density 
perturbations and cosmological parameters. 

We have been able to analyze the effects on the radiation density of the universe of 
new light degrees of freedom which decouple after the QCD phase transition. This includes 
species that decouple at 'complicated' cosmological times, such as the time around which the 
muon becomes non-relativistic. We are able to compute the energy density, and consequently 
Ag*, to an accuracy of 1%. This allows us to place constraints on the couplings in those 
well-motivated BSM effective models which contain new light degrees of freedom, which 
are competitive with constraints coming from other areas of physics. We do this using 
a program which solves the Boltzmann and Friedmann equations for the case of one new 
light species, calculating the resulting evolution of that species' distribution function, while 
approximating the SM species using fully thermalized equilibrium distributions and only 
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Dirac: A > 5 TeV 


5%; see text) 
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Llectroweak Interactions 


T 1 1 1 J_ 

Inconclusive; mass- dependent 


u(iy 


eexAx 


e < 10" 8 for 10 MeV < m x < 150 MeV 
m x > 150 MeV: Decouples during/before 
QCD phase transition 


A'-dipole 




Flavor-blind: Already excluded 
Muon-only: A > 10 3 TeV 



TABLE IV: Compatibility of those natural, minimal models considered here with the recent results 
of the Planck satellite, 5* = 3.50±0.12 and N eff = 3.30±0.27 [21]. While the current Planck results 
are in tension with other observational measurements, future experiments will greatly improve the 
precision and reach of these exclusion limits. 



considering the effects of leading order interaction terms. Using these calculations, we have 
demonstrated the ability of Planck and future experiments to place exclusion limits on all 
natural, minimal models with new light species. The compatibility of each model with the 
recent Planck results is given in Table IIVI 

Higher levels of calculational accuracy could be achieved if we used a different numerical 
algorithm which was better adapted for the integro-differential equations considered in this 
paper, or used a larger and finer momentum grid. In addition, loop corrections to the 
amplitudes, three-body final states, and finite-temperature QFT effects all contribute at the 
0.1% level. If much higher precision is ever achieved observationally, potentially through 
next-generation polarization measurements, then these improvements would be warranted. 
Such a high-precision measurement of g* would better reveal degrees of freedom which 
decouple before or during the QCD phase transition. In such a scenario, this measurement, 
combined with independent measurements of the nature and couplings of a new light degree 
of freedom could potentially even allow us, in this way, to probe the structure of the QCD 
phase transition. 

The future work we intend to pursue is the inclusion of the mass effects on different 
observables in the Cosmic Microwave Background. While this is only relevant in a narrow 
mass range (close to recombination temperatures), it turns out to be quite important for a 
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number of specific models of sterile neutrinos. The more accurately we can describe their 
impact on the ISW effect and on Silk damping, the greater the possibility of finding a 
'smoking gun' for such models. 

If we coarsely divide the types of possible undiscovered particles into four types, catego- 
rized by stable or unstable and light or heavy, this work is an attempt to help push forward 
our probe of one of these types - new stable light particles. As the challenge to build new, 
more powerful high-energy colliders intensifies, it is exciting to see this new frontier mature 
as an additional source of information about the world beyond the Standard Model. 



Appendix A: Notations and conventions 

We take the metric signature to be (+,—,—,—). We set h = c = ks = 1 and give 
temperatures in units of energy. We use /„. ~ 93 MeV. We use M p i « 1.22 x 10 19 GeV, with 
G = M^ 2 . a(t) represents the scale factor of the universe. We use both Weyl and Dirac 
notation for fermions, depending on need. Where there is an ambiguity, we use capital \1/ 
and X for Dirac spinors, and lowercase Greek and x s with subscripts l and C R for Weyl 
spinors. All of our Weyl spinors are left-handed. We use the notation of (the mostly minus 



version of) [62] for two-component spinor Feynman rules. The SM Higgs doublet is called 
h, and SM fermions are generically referred to as ip, when not referencing a specific fermion. 
All BSM light fields not coming from special UV completions are given by for scalars, x 
for fermions and A' for gauge bosons, with associated field strength A'^ v '. 



Appendix B: Details of code 

In order to solve the Boltzmann equation to obtain the distribution functions and con- 
sequently Ap*, a numerical code was written in C++. The code evolves a universe forward 
in time subject to some initial boundary conditions. The inputs to the code are the masses 
uii, statistics er, and degrees of freedom of species g% in the universe, as well as the initial 
temperature Tj. Interactions between the various species in the SM sector are strong enough 
that it is safe to assume that the distribution functions are Fermi-Dirac or Bose-Einstein 
until down to temperatures well below their mass, as discussed in subsection III B[ At that 
point, their number density has become low enough that their interactions to our new species 
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have frozen out, and we work in an effectively radiation-dominated universe, so the error in 
the distribution function does not affect our evolution. In order to work with 0(1) numbers 
for the distribution functions which are decoupling, we track v(p,t) = v instead, defined 



The code solves for the following quantities: 

• Vi(p,t) for all BSM species % 

• Tsm = T~f = T (as we work above the neutrino decoupling temperature) 

• pi, Pi, rii for all species i 

• H and a 

The following equations are used to solve for the aforementioned quantities: 



implicitly through the equation f(p,t) 



1 



where a is 1 for bosons and —1 for fermions. 




df 



dv 



C[f], 



(Bl) 



where we discuss computation of C[f] below, 



dptot 
dt 

where y SM = Y.icSMVi and 



= Vsm 



OTsm Qp x 



dt dt 



3H(p tot +p tot ), 



(B2) 




(B3) 




(B4) 




(B5) 




(B6) 




(B7) 
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2tt 2 



/ dp p 2 f, 



1 1 



(B8) 



30p, 



(B9) 



9*,i 



7T 2 T 4 ' 



The various quantities are tracked over 2000 timesteps, spaced logarithmically. We begin 
at T = 200 MeV and typically end at around 1 MeV. As we only have earlier time infor- 
mation, time derivatives that cannot be computed analytically are typically computed by 
forming an interpolating polynomial to the previous four pieces of data and taking an exact 
derivative of the resulting polynomial. Furthermore, this technique is also used to obtain an 
estimate of the next value of the variable in question, in order to improve the accuracy of the 
code. We found that this technique is considerably more accurate at numerically evaluating 
derivatives than more elementary finite difference methods. 

Our distribution function is evaluated on a grid of 100 momentum-steps, logarithmically 
spaced between 10 keV and 10 GeV. All derivatives with respect to momentum that cannot 
be computed analytically are computed with the interpolating polynomials method described 
previously. 

The algorithm used is as follows: 

• Set boundary conditions: Bose-Einstein or Fermi-Dirac distributions for every species 
with temperature T i7 and a, = 1 

• Compute all initial p iy p iy rij, H, g* ti 

• Compute the next SM temperature with equation flB2j) 

Iterate the following at timestep j until 2000 timesteps: 

• Approximate Hj at the next timestep by extrapolating the previous values of H(t) 

• Solve the Boltzmann equations of all BSM species i for Vij, described more thoroughly 
below 

• Compute all remaining undetermined parameters, as well as H(tj) 

• Compute the next SM temperature with equation (1B2[) 
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The Boltzmann equation (IBlj) was solved using a generalization of a predictor-corrector 
method. The objective was to simultaneously vary at all points on the momentum grid, 
attempting to minimize the quantity 



100 

^-RQ1\/T ; 1 " K 



Ek 



dt 



Hp k E, 



dpk 



pMfiiPk)) 



(BIO) 



iCBSM fc=l 

Because the collisional integral is the most computationally intensive part of the algo- 
rithm, we attempted to minimize the number of calls of it. This was accomplished by 
primarily studying the effects of the variation of t>j on the left-hand side of the Boltzmann 
equation, as the right-hand side varied more slowly, only recomputing C when we had settled 
on a v that minimized the local error 



(Bll) 



dt dp k df t 

with a relative accuracy of ~ 1CT 6 . 

The collisional integral in equation ([6]) has been computed by following the method 



devised by 



We briefly summarize the algorithm; see [8] for more details. Define the 



following quantities: 

• The angle between p\ and p 2 is a 

• The angle between p\ and p 3 is 6 

• The azimuthal angle between p 2 and p 3 is 

• x = cos a 

• z = cos 9 

• Q = ml + m\ + m\ — m\ 

In the amplitude \A4\ 2 , we plug in 



pt -p 2 = EiE 2 - \pi\\p 2 \x, 
pi -p 3 = EiE z - \pi\\p 3 \z, 



p 1 .p 4 = m 2 l + ExE 2 - \pi\\p 2 \x - EiE 3 + \pi\\p 3 \z, 



(B12) 
(B13) 
(B14) 



69 



p 2 -p 3 = E 1 E 2 - \p1Wp2\x - E 1 E 3 + \pi\\p 3 \z + Q/2, (B15) 

p 2 ■ p 4 = E l E 3 - \p x \\p 3 \z + m 2 - Q/2, (B16) 

p 3 ■ p A = E X E 2 - \p x \\p 2 \x - m\ + Q/2. (B17) 

We can change variables to \p\\, \p 2 \, \p 3 \, p~l, x, z, (3 and ji, where n is an integration 
variable parameterizing the S0(2) rotational symmetry about p\. The collisional integral 
has no dependence on /z, and it can therefore be done trivially. After using the momentum- 
conserving delta function to integrate p^, we can use the energy- conserving delta function to 
integrate (3. Now that there are no more four- vectors in our expression, we switch notation 
\Pi\ — > Pi- After some algebra, it has been shown that C can be written in the form 

c(/(pi)) = f d P2 r d P3 / ^ f ji , (bis) 



/o "Jo ^ (27r) 5 16£ 2 £ 3 ' 
where f2 was defined before as f 3 fi(l + 0"i/i)(l + 02/2) — fif 2 (l + cr 3 /3)(l + 04/4) and 
F = F{p u p 2 ,p 3 ) is 



F=jdzj dx Jf'^f^ , MA), (B19) 
y/a(z)x /i + o{z)x + c{z) 



1 r x + 1 n/r\2t 

1 Jx. 



where 



a(z) = {-ApKpI + p\)) + {8p 2 2 p lP3 ) z, (B20) 

b(z) = (pip 2 (8 7 + 4Q)) + p 2 p 3 (8pl - 87 - 4Q) z + {-8 Pl p 2 pl) z 2 , (B21) 
c(z) = (Apjpj - 4 7 2 - 4 7 Q - Q 2 ) + (-Pip 3 (8 7 + 4Q)) 2 + (-4p 2 (p 2 + p 2 )) ^ 2 , (B22) 

- f = E 1 E 2 -E 3 (E 1 + E 2 ), (B23) 

x ± = (B24) 
Note that since a < 0, we know that x + > x_. If we define 

1 



z± = 



2pm (~ 27 - 2p i- Q± 2p2 \I 2 ^ + Pi+Pl+Pl + Q) > (B25) 
then 0(A) is 1 when 2 7 + + p\ + p| + Q > 0, z + > —1, and z- < 1, and is otherwise. 
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The amplitudes \Ai\ 2 were computed with the assistance of Tracer 108]. After the 
substitutions above, \Ai\ 2 can be written as a rational function in x and z. We first expand 
it in x, and integrate it analytically with Mathematica. Afterwards, if it is possible to 
analytically integrate with respect to z, we do so and store the results at all 10 6 combinations 
{Pi,P2,Ps}- Otherwise, we numerically integrate with respect to z at all 10 6 points in the 
resulting phase space. These are stored and then loaded into our C++ code. 



The code has been verified to give the same answer for T e jf :l/ (p) as that given in jfj|, giving 
the same values for Ag* to the percent level. In addition, the results were computed and 
compared to all cases where it is possible to use the instantaneous decoupling approximation 
or otherwise solve the problem analytically, and agreement was again found to the percent 
level or better in all cases. Percent-level accuracy is more precise than the resolution of the 
Planck satellite, and so we do not quote theoretical errors throughout the paper. 

In subsection IIII C 1\ we reference a modified version of the code suitable for tracking 
two separate thermalized sectors which undergo partially thermalizing interactions. The 
structure of the code is very similar to the code outlined above, but with a few changes: 

• The initial conditions are Tsm, Thid and e. 

• The distribution functions are always Bose-Einstein or Fermi-Dirac, and so only the 
two temperatures are tracked, eliminating the need for storing any distribution func- 
tions. 

• We use a modified set of equations to track the distribution functions: 



dT SM -3H(p SM + Psm) ~ T , . 

at y SM 

dT hid -3H(p hi d + Phid) + r 



dt y hid 



(B27) 



9 <x 



2tt 2 



/■oo 

/ dp x p\ C[h{ Pl )\. (B28) 
Jo 



A semi-implicit Euler method was used to compute the temperatures at the next 
timesteps, in order to minimize the amount of time spent computing V. 
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